home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / test < prev    next >
Encoding:
Text File  |  1994-11-09  |  68.6 KB  |  3,089 lines  |  [TEXT/ttxt]
  1. tic();
  2. //
  3. //  RLaB test program
  4. //
  5.  
  6. //
  7. // Test Parameters and some functions we will need later
  8. //
  9.  
  10. rand("default");
  11. pi = 4.0*atan(1.0);
  12. X = 3;
  13.  
  14. epsilon = function() 
  15. {
  16.   eps = 1.0;
  17.   while((eps + 1.0) != 1.0) 
  18.   {
  19.     eps = eps/2.0;
  20.   }
  21.   return eps;
  22. };
  23.  
  24. eps = epsilon();
  25.  
  26. eye = function( m , n ) 
  27. {
  28.  
  29.   if (!exist (n))
  30.   {
  31.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  32.     new = zeros (m[1], m[2]);
  33.     N = min ([m[1], m[2]]);
  34.   else
  35.     if (class (m) == "string" || class (n) == "string") {
  36.       error ("eye(), string arguments not allowed");
  37.     }
  38.     if (max (size (m)) == 1 && max (size (n)) == 1)
  39.     {
  40.       new = zeros (m[1], n[1]);
  41.       N = min ([m[1], n[1]]);
  42.     else
  43.       error ("matrix arguments to eye() must be 1x1");
  44.     }
  45.   }
  46.   for(i in 1:N)
  47.   {
  48.     new[i;i] = 1.0;
  49.   }
  50.   return new;
  51. };
  52.  
  53. symm = function( A )
  54. {
  55.   return (A + A')./2;
  56. };
  57.  
  58. //-------------------------------------------------------------------//
  59.  
  60. // Synopsis:    Pascal matrix.
  61.  
  62. // Syntax:  P = pascal ( N )
  63.  
  64. // Description:
  65.  
  66. //  The Pascal matrix of order N: a symmetric positive definite
  67. //  matrix with integer entries taken from Pascal's triangle. The
  68. //  Pascal matrix is totally positive and its inverse has integer
  69. //  entries.  Its eigenvalues occur in reciprocal pairs. COND(P)
  70. //  is approximately 16^N/(N*PI) for large N. PASCAL(N,1) is the
  71. //  lower triangular Cholesky factor (up to signs of columns) of
  72. //  the Pascal matrix.   It is involutary (is its own
  73. //  inverse). PASCAL(N,2) is a transposed and permuted version of
  74. //  PASCAL(N,1) which is a cube root of the identity.
  75.  
  76. //      References:
  77. //      R. Brawer and M. Pirovino, The linear algebra of the Pascal matrix,
  78. //           Linear Algebra and Appl., 174 (1992), pp. 13-23 (this paper
  79. //           gives a factorization of L = PASCAL(N,1) and a formula for the
  80. //           elements of L^k).
  81. //      S. Karlin, Total Positivity, Volume 1, Stanford University Press,
  82. //           1968.  (Page 137: shows i+j-1 choose j is TP (i,j=0,1,...).
  83. //                   PASCAL(N) is a submatrix of this matrix.)
  84. //      M. Newman and J. Todd, The evaluation of matrix inversion programs,
  85. //           J. Soc. Indust. Appl. Math., 6(4):466--476, 1958.
  86. //      H. Rutishauser, On test matrices, Programmation en Mathematiques
  87. //           Numeriques, Editions Centre Nat. Recherche Sci., Paris, 165,
  88. //           1966, pp. 349-365.  (Gives an integral formula for the
  89. //           elements of PASCAL(N).)
  90. //      J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  91. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 172.
  92. //      H.W. Turnbull, The Theory of Determinants, Matrices, and Invariants,
  93. //           Blackie, London and Glasgow, 1929.  (PASCAL(N,2) on page 332.)
  94.  
  95. //  This file is a translation of pascal.m from version 2.0 of
  96. //  "The Test Matrix Toolbox for Matlab", described in Numerical
  97. //  Analysis Report No. 237, December 1993, by N. J. Higham.
  98.  
  99. //-------------------------------------------------------------------//
  100.  
  101. pascal = function ( n , k )
  102. {
  103.   local(P, i, j, k, n)
  104.  
  105.   if (!exist (k)) { k = 0; }
  106.  
  107.   P = diag( (-1).^[0:n-1] );
  108.   P[; 1] = ones(n,1);
  109.  
  110.   //  Generate the Pascal Cholesky factor (up to signs).
  111.  
  112.   for (j in 2:n-1)
  113.   {
  114.     for (i in j+1:n)
  115.     {
  116.       P[i;j] = P[i-1;j] - P[i-1;j-1];
  117.     }
  118.   }
  119.  
  120.   if (k == 0)
  121.   {
  122.     P = P*P';
  123.   else if (k == 2) {
  124.     P = rot90(P,3);
  125.     if (n/2 == round(n/2)) { P = -P; }
  126.   }}
  127.  
  128.   return P;
  129. };
  130.  
  131. //-------------------------------------------------------------------//
  132.  
  133. // Synopsis:    Random, orthogonal upper Hessenberg matrix.
  134.  
  135. // Syntax:      H = ohess ( N )
  136.  
  137. // Description:
  138.  
  139. //      H is an N-by-N real, random, orthogonal upper Hessenberg
  140. //      matrix. Alternatively, H = OHESS(X), where X is an arbitrary
  141. //      real N-vector (N > 1) constructs H non-randomly using the
  142. //      elements of X as parameters. In both cases H is constructed
  143. //      via a product of N-1 Givens rotations. 
  144.  
  145. //      Note: See Gragg (1986) for how to represent an N-by-N
  146. //      (complex) unitary Hessenberg matrix with positive subdiagonal
  147. //      elements in terms of 2N-1 real parameters (the Schur
  148. //      parametrization). This M-file handles the real case only and
  149. //      is intended simply as a convenient way to generate random or
  150. //      non-random orthogonal Hessenberg matrices.
  151.  
  152. //      Reference:
  153. //      W.B. Gragg, The QR algorithm for unitary Hessenberg matrices,
  154. //      J. Comp. Appl. Math., 16 (1986), pp. 1-8.
  155.  
  156. //  This file is a translation of ohess.m from version 2.0 of
  157. //  "The Test Matrix Toolbox for Matlab", described in Numerical
  158. //  Analysis Report No. 237, December 1993, by N. J. Higham.
  159.  
  160. //-------------------------------------------------------------------//
  161.  
  162. ohess = function ( x )
  163. {
  164.   local (x, n, H, theta, c, s, i)
  165.   global (pi)
  166.  
  167.   if (any (imag (x))) { error("Parameter must be real."); }
  168.  
  169.   n = max(size(x));
  170.  
  171.   if (n == 1)
  172.   {
  173.     //  Handle scalar x.
  174.     n = x;
  175.     x = rand(n-1, 1)*2*pi;
  176.     H = eye(n,n);
  177.     H[n;n] = sign(rand());
  178.   else
  179.     H = eye(n,n);
  180.     H[n;n] = sign(x[n]) + (x[n]==0);   // Second term ensures H[n;n] nonzero.
  181.   }
  182.  
  183.   for (i in n:2:-1)
  184.   {
  185.     // Apply Givens rotation through angle x(i-1).
  186.     theta = x[i-1];
  187.     c = cos(theta);
  188.     s = sin(theta);
  189.     H[ [i-1, i] ;] = [ c*H[i-1;]+s*H[i;] ;
  190.                        -s*H[i-1;]+c*H[i;] ];
  191.   }
  192.  
  193.   return H;
  194. };
  195.  
  196. //-------------------------------------------------------------------//
  197.  
  198. // Synopsis:    Random matrix with pre-assigned singular values.
  199.  
  200. // Syntax:  R = randsvd (N, KAPPA, MODE, KL, KU)
  201.  
  202. // Description:
  203.  
  204. //  R is a (banded) random matrix of order N with COND(A) = KAPPA
  205. //  and singular values from the distribution MODE.
  206.  
  207. //      N may be a 2-vector, in which case the matrix is N(1)-by-N(2).
  208. //      Available types:
  209. //             MODE = 1:   one large singular value,
  210. //             MODE = 2:   one small singular value,
  211. //             MODE = 3:   geometrically distributed singular values,
  212. //             MODE = 4:   arithmetically distributed singular values,
  213. //             MODE = 5:   random singular values with unif. dist. logarithm.
  214.  
  215. //      If omitted, MODE defaults to 3, and KAPPA defaults to SQRT(1/EPS).
  216. //      If MODE < 0 then the effect is as for ABS(MODE) except that in the
  217. //      original matrix of singular values the order of the diagonal entries
  218. //      is reversed: small to large instead of large to small.
  219. //      KL and KU are the lower and upper bandwidths respectively; if they
  220. //      are omitted a full matrix is produced.
  221. //      If only KL is present, KU defaults to KL.
  222. //      Special case: if KAPPA < 0 then a random full symmetric positive
  223. //                    definite matrix is produced with COND(A) = -KAPPA and
  224. //                    eigenvalues distributed according to MODE.
  225. //                    KL and KU, if present, are ignored.
  226.  
  227. //      This routine is similar to the more comprehensive Fortran routine xLATMS
  228. //      in the following reference:
  229. //      J.W. Demmel and A. McKenney, A test matrix generation suite,
  230. //      LAPACK Working Note #9, Courant Institute of Mathematical Sciences,
  231. //      New York, 1989.
  232.  
  233. //  This file is a translation of randsvd.m from version 2.0 of
  234. //  "The Test Matrix Toolbox for Matlab", described in Numerical
  235. //  Analysis Report No. 237, December 1993, by N. J. Higham.
  236.  
  237. // Dependencies
  238. #   rfile bandred
  239. #   rfile qmult
  240.  
  241. //-------------------------------------------------------------------//
  242.  
  243. randsvd = function (n, kappa, mode, kl, ku)
  244. {
  245.   local (n, kappa, mode, kl, ku, A, sigma, ...
  246.          posdef, p, m, j, Q, factor)
  247.   global(qmult)
  248.  
  249.   if (!exist (kappa)) { kappa = sqrt(1/eps); }
  250.   if (!exist (mode))  { mode = 3; }
  251.   if (!exist (kl))    { kl = n-1; } // Full matrix.
  252.   if (!exist (ku))    { ku = kl; }  // Same upper and lower bandwidths.
  253.  
  254.   if (abs(kappa) < 1) {
  255.     error("Condition number must be at least 1!");
  256.   }
  257.  
  258.   posdef = 0;
  259.   if (kappa < 0)            // Special case.
  260.   {
  261.     posdef = 1;
  262.     kappa = -kappa;
  263.   }
  264.  
  265.   p = min(n);
  266.   m = n[1];     // Parameter n specifies dimension: m-by-n.
  267.   n = n[max(size(n))];
  268.  
  269.   if (p == 1)       // Handle case where A is a vector.
  270.   {
  271.     rand("normal", 0, 1);
  272.     A = rand(m, n);
  273.     A = A/norm(A, "2");
  274.     return A;
  275.   }
  276.  
  277.   j = abs(mode);
  278.  
  279.   // Set up vector sigma of singular values.
  280.   if (j == 3)
  281.   {
  282.     factor = kappa^(-1/(p-1));
  283.     sigma = factor.^[0:p-1];
  284.  
  285.   else if (j == 4) {
  286.     sigma = ones(p,1) - (0:p-1)'/(p-1)*(1-1/kappa);
  287.  
  288.   else if (j == 5) {    // In this case cond(A) <= kappa.
  289.     rand("uniform", 0, 1)
  290.     sigma = exp( -rand(p,1)*log(kappa) );
  291.  
  292.   else if (j == 2) {
  293.     sigma = ones(p,1);
  294.     sigma[p] = 1/kappa;
  295.  
  296.   else if (j == 1) {
  297.     sigma = ones(p,1)./kappa;
  298.     sigma[1] = 1;
  299.  
  300.   }}}}}
  301.  
  302.  
  303.   // Convert to diagonal matrix of singular values.
  304.   if (mode < 0) {
  305.     sigma = sigma[p:1:-1];
  306.   }
  307.  
  308.   sigma = diag(sigma);
  309.  
  310.   if (posdef)       // Handle special case.
  311.   {
  312.     Q = qmult(p);
  313.     A = Q'*sigma*Q;
  314.     A = (A + A')/2; // Ensure matrix is symmetric.
  315.     return A;
  316.   }
  317.  
  318.   if (m != n) 
  319.   {
  320.     sigma[m; n] = 0;    // Expand to m-by-n diagonal matrix.
  321.   }
  322.  
  323.   if (kl == 0 && ku == 0)   // Diagonal matrix requested - nothing more to do.
  324.   {
  325.     A = sigma;
  326.     return A;
  327.   }
  328.  
  329.   // A = U*sigma*V, where U, V are random orthogonal matrices from the
  330.   // Haar distribution.
  331.   A = qmult(sigma');
  332.   A = qmult(A');
  333.  
  334.   if (kl < n-1 || ku < n-1) // Bandwidth reduction.
  335.   {
  336.    A = bandred(A, kl, ku);
  337.   }
  338.  
  339.   rand("default");
  340.   return A;
  341. };
  342.  
  343. //-------------------------------------------------------------------//
  344.  
  345. // Synopsis:    Band reduction by two-sided unitary transformations.
  346.  
  347. // Syntax:  bandred ( A , KL, KU )
  348.  
  349. // Description:
  350.  
  351. //  bandred(A, KL, KU) is a matrix unitarily equivalent to A with
  352. //  lower bandwidth KL and upper bandwidth KU (i.e. B(i,j) = 0 if
  353. //  i > j+KL or j > i+KU).  The reduction is performed using
  354. //  Householder transformations. If KU is omitted it defaults to
  355. //  KL. 
  356.  
  357. //  Called by randsvd.
  358. //  This is a `standard' reduction.  Cf. reduction to bidiagonal
  359. //  form prior to computing the SVD.  This code is a little
  360. //  wasteful in that it computes certain elements which are
  361. //  immediately set to zero! 
  362. //
  363. //      Reference:
  364. //      G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  365. //      Johns Hopkins University Press, Baltimore, Maryland, 1989.
  366. //      Section 5.4.3.
  367.  
  368. //  This file is a translation of bandred.m from version 2.0 of
  369. //  "The Test Matrix Toolbox for Matlab", described in Numerical
  370. //  Analysis Report No. 237, December 1993, by N. J. Higham.
  371.  
  372. // Dependencies
  373. #   rfile house
  374.  
  375. //-------------------------------------------------------------------//
  376.  
  377. bandred = function ( A , kl , ku )
  378. {
  379.   local (A, kl, ku, beta, v, flip, ltmp, j, temp, m, n)
  380.  
  381.   if (!exist (ku)) { ku = kl; else ku = ku; }
  382.  
  383.   if (kl == 0 && ku == 0) {
  384.     error ("You''ve asked for a diagonal matrix.  In that case use the SVD!");
  385.   }
  386.  
  387.   // Check for special case where order of left/right transformations matters.
  388.   // Easiest approach is to work on the transpose, flipping back at the end.
  389.  
  390.   flip = 0;
  391.   if (ku == 0)
  392.   {
  393.     A = A';
  394.     temp = kl; kl = ku; ku = temp; flip = 1;
  395.   }
  396.  
  397.   m = A.nr; n = A.nc; 
  398.  
  399.   for (j in 1 : min( min(m, n), max(m-kl-1, n-ku-1) ))
  400.   {
  401.     if (j+kl+1 <= m)
  402.     {
  403.        ltmp = house(A[j+kl:m;j]);
  404.        beta = ltmp.beta; v = ltmp.v;
  405.        temp = A[j+kl:m;j:n];
  406.        A[j+kl:m;j:n] = temp - beta*v*(v'*temp);
  407.        A[j+kl+1:m;j] = zeros(m-j-kl,1);
  408.     }
  409.  
  410.     if (j+ku+1 <= n)
  411.     {
  412.        ltmp = house(A[j;j+ku:n]');
  413.        beta = ltmp.beta; v = ltmp.v;
  414.        temp = A[j:m;j+ku:n];
  415.        A[j:m;j+ku:n] = temp - beta*(temp*v)*v';
  416.        A[j;j+ku+1:n] = zeros(1,n-j-ku);
  417.     }
  418.   }
  419.  
  420.   if (flip) {
  421.     A = A';
  422.   }
  423.  
  424.   return A;
  425. };
  426.  
  427. //-------------------------------------------------------------------//
  428.  
  429. // Synopsis:    Lehmer matrix - symmetric positive definite.
  430.  
  431. // Syntax:      A = lehmer ( N )
  432.  
  433. // Description:
  434.  
  435. //      A is the symmetric positive definite N-by-N matrix with
  436. //                     A(i,j) = i/j for j >= i.
  437. //      A is totally nonnegative.  INV(A) is tridiagonal, and explicit
  438. //      formulas are known for its entries. 
  439.  
  440. //      N <= COND(A) <= 4*N*N.
  441.  
  442. //      References:
  443. //        M. Newman and J. Todd, The evaluation of matrix inversion
  444. //           programs, J. Soc. Indust. Appl. Math., 6 (1958), pp. 466-476.
  445. //        Solutions to problem E710 (proposed by D.H. Lehmer): The inverse
  446. //           of a matrix, Amer. Math. Monthly, 53 (1946), pp. 534-535.
  447. //        J. Todd, Basic Numerical Mathematics, Vol. 2: Numerical Algebra,
  448. //           Birkhauser, Basel, and Academic Press, New York, 1977, p. 154.
  449.  
  450. //  This file is a translation of lehmer.m from version 2.0 of
  451. //  "The Test Matrix Toolbox for Matlab", described in Numerical
  452. //  Analysis Report No. 237, December 1993, by N. J. Higham.
  453.  
  454. //-------------------------------------------------------------------//
  455.  
  456. lehmer = function ( n )
  457. {
  458.   local (A)
  459.   global (tril)
  460.  
  461.   A = ones(n,1)*(1:n);
  462.   A = A./A';
  463.   A = tril(A) + tril(A,-1)';
  464.  
  465.   return A;
  466. };
  467.  
  468. //-------------------------------------------------------------------//
  469.  
  470. // Synopsis:    Pre-multiply by random orthogonal matrix.
  471.  
  472. // Syntax:  B = qmult ( A )
  473.  
  474. // Description:
  475.  
  476. //  B is Q*A where Q is a random real orthogonal matrix from the
  477. //  Haar distribution, of dimension the number of rows in
  478. //  A. Special case: if A is a scalar then QMULT(A) is the same as
  479.  
  480. //      qmult(eye(a))
  481.  
  482. //       Called by RANDSVD.
  483.  
  484. //       Reference:
  485. //       G.W. Stewart, The efficient generation of random
  486. //       orthogonal matrices with an application to condition estimators,
  487. //       SIAM J. Numer. Anal., 17 (1980), 403-409.
  488.  
  489. //  This file is a translation of qmult.m from version 2.0 of
  490. //  "The Test Matrix Toolbox for Matlab", described in Numerical
  491. //  Analysis Report No. 237, December 1993, by N. J. Higham.
  492.  
  493. //-------------------------------------------------------------------//
  494.  
  495. qmult = function ( A )
  496. {
  497.   local (A, B, n, m, d, k, x, s, sgn, y, i, beta)
  498.  
  499.   n = A.nr; m = A.nc;
  500.  
  501.   //  Handle scalar A.
  502.   if (max(n,m) == 1)
  503.   {
  504.     n = A;
  505.     A = eye(n,n);
  506.   }
  507.  
  508.   d = zeros(n,n);
  509.  
  510.   for (k in n-1:1:-1)
  511.   {
  512.     // Generate random Householder transformation.
  513.     rand("normal", 0, 1);
  514.     x = rand(n-k+1,1);
  515.     s = norm(x, "2");
  516.     sgn = sign(x[1]) + (x[1]==0);   // Modification for sign(1)=1.
  517.     s = sgn*s;
  518.     d[k] = -sgn;
  519.     x[1] = x[1] + s;
  520.     beta = s*x[1];
  521.  
  522.     // Apply the transformation to A.
  523.     y = x'*A[k:n;];
  524.     A[k:n;] = A[k:n;] - x*(y/beta);
  525.   }
  526.  
  527.   // Tidy up signs.
  528.   for (i in 1:n-1)
  529.   {
  530.     A[i;] = d[i]*A[i;];
  531.   }
  532.  
  533.   A[n;] = A[n;]*sign(rand());
  534.   B = A;
  535.  
  536.   rand("default");
  537.   return B;
  538. };
  539.  
  540. //-------------------------------------------------------------------//
  541.  
  542. // Synopsis:    Create a Householder matrix.
  543.  
  544. // Syntax:  house ( X )
  545.  
  546. //  If HOUSE(x), which returns a list containing elements `v' and
  547. //  `beta',  then H = EYE - beta*v*v' is a Householder matrix such
  548. //  that Hx = -sign(x(1))*norm(x)*e_1. 
  549. //  NB: If x = 0 then v = 0, beta = 1 is returned.
  550. //            x can be real or complex.
  551. //            sign(x) := exp(i*arg(x)) ( = x./abs(x) when x ~= 0).
  552.  
  553. //  Theory: (textbook references Golub & Van Loan 1989, 38-43;
  554. //           Stewart 1973, 231-234, 262; Wilkinson 1965, 48-50).
  555. //      Hx = y: (I - beta*v*v')x = -s*e_1.
  556. //      Must have |s| = norm(x), v = x+s*e_1, and
  557. //      x'y = x'Hx =(x'Hx)' real => arg(s) = arg(x(1)).
  558. //      So take s = sign(x(1))*norm(x) (which avoids cancellation).
  559. //      v'v = (x(1)+s)^2 + x(2)^2 + ... + x(n)^2
  560. //          = 2*norm(x)*(norm(x) + |x(1)|).
  561.  
  562. //  References:
  563. //        G.H. Golub and C.F. Van Loan, Matrix Computations, second edition,
  564. //           Johns Hopkins University Press, Baltimore, Maryland, 1989.
  565. //        G.W. Stewart, Introduction to Matrix Computations, Academic Press,
  566. //           New York, 1973,
  567. //        J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University
  568. //           Press, 1965.
  569.  
  570. //  This file is a translation of house.m from version 2.0 of
  571. //  "The Test Matrix Toolbox for Matlab", described in Numerical
  572. //  Analysis Report No. 237, December 1993, by N. J. Higham.
  573.  
  574. //-------------------------------------------------------------------//
  575.  
  576. house = function ( x )
  577. {
  578.   local (m, n, s, v, beta)
  579.  
  580.   m = x.nr; n = x.nc;
  581.   if (n > 1) { error ("Argument must be a column vector."); }
  582.  
  583.   s = norm(x,"2") * (sign(x[1]) + (x[1]==0)); // Modification for sign(0)=1.
  584.   v = x;
  585.   if (s == 0)   // Quit if x is the zero vector.
  586.   {
  587.     beta = 1; 
  588.     return << beta = beta ; v = v >>;
  589.   }
  590.  
  591.   v[1] = v[1] + s;
  592.   beta = 1/(s'*v[1]);       // NB the conjugated s.
  593.  
  594.   // beta = 1/(abs(s)*(abs(s)+abs(x(1)) would guarantee beta real.
  595.   // But beta as above can be non-real (due to rounding) only 
  596.   // when x is complex.
  597.  
  598.   return << beta = beta ; v = v >>;
  599. };
  600.  
  601. //-------------------------------------------------------------------//
  602.  
  603. //  Syntax: tril ( A )
  604. //      tril ( A , K )
  605.  
  606. //  Description:
  607.  
  608. //  tril(x) returns the lower triangular part of A.
  609.  
  610. //  tril(A,K) returns the elements on and below the K-th diagonal of
  611. //  A.
  612.  
  613. //  K = 0: main diagonal
  614. //  K > 0: above the main diag.
  615. //  K < 0: below the main diag.
  616.  
  617. //  See Also: triu
  618. //-------------------------------------------------------------------//
  619.  
  620. tril = function(x, k) 
  621. {
  622.   local(i, j, nr, nc, y);
  623.  
  624.   if (!exist (k)) { k = 0; }
  625.   nr = x.nr; nc = x.nc;
  626.   if(k > 0) 
  627.   { 
  628.     if (k > (nc - 1)) { error ("tril: invalid value for k"); }
  629.   else
  630.     if (abs (k) > (nr - 1)) { error ("tril: invalid value for k"); }
  631.   }
  632.  
  633.   y = zeros(nr, nc);
  634.  
  635.   for(i in max( [1,1-k] ):nr) {
  636.     j = 1:min( [nc, i+k] );
  637.     y[i;j] = x[i;j];
  638.   }
  639.  
  640.   return y;
  641. };
  642.  
  643. //-------------------------------------------------------------------//
  644.  
  645. //  Syntax: triu ( A )
  646. //      triu ( A , K )
  647.  
  648. //  Description:
  649.  
  650. //  triu(x) returns the upper triangular part of A.
  651.  
  652. //  tril(x; k) returns the elements on and above the k-th diagonal of
  653. //  A. 
  654.  
  655. //  K = 0: main diagonal
  656. //  K > 0: above the main diag.
  657. //  K < 0: below the main diag.
  658.  
  659. //  See Also: tril
  660. //-------------------------------------------------------------------//
  661.  
  662. triu = function(x, k) 
  663. {
  664.   local(i, j, nr, nc, y);
  665.  
  666.   if (!exist (k)) { k = 0; }
  667.   nr = x.nr; nc = x.nc;
  668.  
  669.   if(k > 0) 
  670.   { 
  671.     if (k > (nc - 1)) { error ("triu: invalid value for k"); }
  672.   else
  673.     if (abs (k) > (nr - 1)) { error ("triu: invalid value for k"); }
  674.   }
  675.  
  676.   y = zeros(nr, nc);
  677.  
  678.   for(j in max( [1,1+k] ):nc) {
  679.     i = 1:min( [nr, j-k] );
  680.     y[i;j] = x[i;j];
  681.   }
  682.  
  683.   return y;
  684. };
  685.  
  686. //
  687. //-------------------- Test Relational Expressions -------------------
  688. //
  689. printf("\tstart scalar tests...\n");
  690. printf("\tstart relational tests...\n");
  691.  
  692. //    SCALAR CONSTANTS (REAL)
  693. if( !(1<2) ) { error(); }
  694. if( !(1<=2) ) { error(); }
  695. if( 1>2 ) { error(); }
  696. if( 1>=2 ) { error(); }
  697. if( 1==2 ) { error(); }
  698. if( !(1!=2) ) { error(); }
  699. if( !1 ) { error(); }
  700. if( !!!1) { error(); }
  701.  
  702. if( !([1]<[2]) ) { error(); }
  703. if( !([1]<=[2]) ) { error(); }
  704. if( [1]>[2] ) { error(); }
  705. if( [1]>=[2] ) { error(); }
  706. if( [1]==[2] ) { error(); }
  707. if( !([1]!=[2]) ) { error(); }
  708. if( ![1] ) { error(); }
  709. if( !!![1]) { error(); }
  710.  
  711. //    SCALAR CONSTANTS (COMPLEX)
  712. if( !(1+2i<2+3i) ) { error(); }
  713. if( !(1+2i<=2+3i) ) { error(); }
  714. if( 1+2i>2+3i ) { error(); }
  715. if( 1+2i>=2+3i ) { error(); }
  716. if( 1+2i==2+3i ) { error(); }
  717. if( !(1+2i!=2+3i) ) { error(); }
  718. if( !1+2i ) { error(); }
  719. if( !!!1+2i) { error(); }
  720.  
  721. if( !([1+2i]<[2+3i]) ) { error(); }
  722. if( !([1+2i]<=[2+3i]) ) { error(); }
  723. if( [1+2i]>[2+3i] ) { error(); }
  724. if( [1+2i]>=[2+3i] ) { error(); }
  725. if( [1+2i]==[2+3i] ) { error(); }
  726. if( !([1+2i]!=[2+3i]) ) { error(); }
  727. if( ![1+2i] ) { error(); }
  728. if( !!![1+2i]) { error(); }
  729.  
  730. //    SCALAR ENTITIES (REAL)
  731. a=1;b=2;
  732. if( !(a<b) ) { error(); }
  733. if( !(a<=b) ) { error(); }
  734. if( a>b ) { error(); }
  735. if( a>=b ) { error(); }
  736. if( a==b ) { error(); }
  737. if( !(a!=b) ) { error(); }
  738. if( !a ) { error(); }
  739. if( !!!a) { error(); }
  740.  
  741. if( !([a]<[b]) ) { error(); }
  742. if( !([a]<=[b]) ) { error(); }
  743. if( [a]>[b] ) { error(); }
  744. if( [a]>=[b] ) { error(); }
  745. if( [a]==[b] ) { error(); }
  746. if( !([a]!=[b]) ) { error(); }
  747. if( ![a] ) { error(); }
  748. if( !!![a]) { error(); }
  749.  
  750. //    SCALAR ENTITIES (COMPLEX)
  751. a=1+2i;b=2+3i;
  752. if( !(a<b) ) { error(); }
  753. if( !(a<=b) ) { error(); }
  754. if( a>b ) { error(); }
  755. if( a>=b ) { error(); }
  756. if( a==b ) { error(); }
  757. if( !(a!=b) ) { error(); }
  758. if( !a ) { error(); }
  759. if( !!!a) { error(); }
  760.  
  761. if( !([a]<[b]) ) { error(); }
  762. if( !([a]<=[b]) ) { error(); }
  763. if( [a]>[b] ) { error(); }
  764. if( [a]>=[b] ) { error(); }
  765. if( [a]==[b] ) { error(); }
  766. if( !([a]!=[b]) ) { error(); }
  767. if( ![a] ) { error(); }
  768. if( !!![a]) { error(); }
  769.  
  770. if (! all (all (-2 <  rand(4,4)))) { error (); }
  771. if (! all (all (-2 <= rand(4,4)))) { error (); }
  772. if (! all (all ( 2 >  rand(4,4)))) { error (); }
  773. if (! all (all ( 2 >= rand(4,4)))) { error (); }
  774.  
  775. if (! all (all (rand(4,4) >  -2))) { error (); }
  776. if (! all (all (rand(4,4) >= -2))) { error (); }
  777. if (! all (all (rand(4,4) <   2))) { error (); }
  778. if (! all (all (rand(4,4) <=  2))) { error (); }
  779.  
  780. if (! all (all (rand (4,4) >  -rand (4,4)))) { error (); }
  781. if (! all (all (rand (4,4) >= -rand (4,4)))) { error (); }
  782. if (! all (all (-rand (4,4) <  rand (4,4)))) { error (); }
  783. if (! all (all (-rand (4,4) <= rand (4,4)))) { error (); }
  784.  
  785. //------------------------- Test REAL SCALARS ------------------------
  786. //
  787. //    CONSTANTS
  788. //      Addition
  789. if(1+2 != 3) { error(); }
  790. //      Subtraction
  791. if(1-2 != -1) { error(); }
  792. //      Multiply
  793. if(1*2 != 2) { error(); }
  794. //      Divide
  795. if(1/2 != 0.5) { error(); }
  796. //      Power
  797. if(2^2 != 4) { error(); }
  798. if(4^0 != 1) { error(); }
  799. //      Unary Minus
  800. if(2-3 != -1) { error(); }
  801. //
  802. //    ENTITIES
  803. //
  804. a = 1; b = 2; c = 3; d = 0.5;
  805. //      Addition
  806. if(a+b != c) { error(); }
  807. //      Subtraction
  808. if(a-b != -a) { error(); }
  809. //      Multiply
  810. if(a*b != b) { error(); }
  811. //      Divide
  812. if(a/b != d) { error(); }
  813. //      Power
  814. if(b^b != b*b) { error(); }
  815. if(b^0 != 1) { error (); }
  816. //      Unary Minus
  817. if(b-c != -a) { error(); }
  818. //
  819. //  ENTITIES & CONSTANTS
  820. //
  821. if(a+2 != c) { error(); }
  822. //      Subtraction
  823. if(1-b != -a) { error(); }
  824. //      Multiply
  825. if(a*2 != 2) { error(); }
  826. //      Divide
  827. if(1/b != d) { error(); }
  828. //      Power
  829. if(2^b != b*b) { error(); }
  830. //      Unary Minus
  831. if(b-3 != -a) { error(); }
  832. //
  833. //------------------------Test COMPLEX SCALARS -------------------------
  834. //
  835. //    CONSTANTS
  836. if(sqrt(-1) != 1i) { error(); }
  837. //      Addition
  838. if((1+2i)+(2+3i) != (3+5i)) { error(); }
  839. //      Subtraction
  840. if((1+2i)-(3+4i) != (-2-2i)) { error(); }
  841. //      Multiply
  842. if((1+2i)*(3+4i) != (-5+10i)) { error(); }
  843. //      Divide
  844. if((1+2i)/(3-4i) != (-.2+.4i)) { error(); }
  845. //      Power
  846. //      Precision problems prevent us from testing these. Have to
  847. //      be checked by hand.
  848. //  (1+2i)^2 = -3 + 4i
  849. //  (1+2i)^.5 = 1.272 + 7.862e-1i
  850. //  if((1+2i)^2 != (-3+4i)) { error(); }
  851. //  if((1+2i)^10 != (237+3116i)) { error(); }
  852. //      Unary Minus
  853. if(-(1+2i) != -1-2i) { error(); }
  854. //
  855. //    ENTITIES
  856. //
  857. a = 1+2i; b = 3+4i; c = 4+6i;
  858. if(a+b != c) { error(); }
  859. //      Subtraction
  860. if(a-b != -2-2i) { error(); }
  861. //      Multiply
  862. if(a*b != -5+10i) { error(); }
  863. //      Divide
  864. //if(a/(3-4i) != -.2+.4i) { error(); }
  865. //      Power
  866. //  if(b^b != b*b) { error(); }
  867. //      Unary Minus
  868. if(-a != -1-2i) { error(); }
  869. //
  870. //    ENTITIES & CONSTANTS
  871. //
  872. if(a+(3+4i) != c) { error(); }
  873. //      Subtraction
  874. if((1+2i)-b != -2-2i) { error(); }
  875. //      Multiply
  876. if(a*(3+4i) != -5+10i) { error(); }
  877. //      Divide
  878. //if(a/(3-4i) != -.2+.4i) { error(); }
  879. //      Power
  880. //if(b^b != b*b) { error(); }
  881. //      Unary Minus
  882. if(-(1+2i) != -1-2i) { error(); }
  883. //
  884. // String - Numerical Equalities
  885. //
  886. if ((1 == "1")) { error(); }
  887. if (([1] == "1")) { error(); }
  888. if (("1" == 1)) { error(); }
  889. if (("1" == [1])) { error(); }
  890.  
  891. if (rand(3,3) == "str") { error(); }
  892. if ("str" == rand(3,3)) { error(); }
  893. if (!any (any ((["1", "2"; "3", "4"] == "1") == [1,0;0,0]))) { error(); }
  894.  
  895. //
  896. //----------------------- Test REAL MATRICES ---------------------------
  897. //
  898.  
  899. printf("\tstart matrix tests...\n");
  900. printf("\t\treal-matrices\n");
  901.  
  902. //  Read in test matrices
  903. //
  904.  
  905. read("test.input");
  906.  
  907. //
  908. //  Matrix construction
  909. //
  910.  
  911. if(any([1;2;3] != [1,2,3]')) {
  912.   error();
  913. }
  914.  
  915. if(any (any (m0 != zeros(2,2)))) { error(); }
  916. if(any (any (m1 != 1+zeros(2,2)))) { error(); }
  917. if(any (any (m2 != [1,2;3,4]))) { error(); }
  918. if(any (any (m3 != [1+2i,2+3i;3+4i,5+6i]))) { error(); }
  919. if(any (any ([1,2;3+0i,4+0i] != m2))) { error(); }
  920. if(any (any ([m2] != m2))) { error(); }
  921.  
  922. //
  923. //  Matrix indexing
  924. //
  925.  
  926. p = pascal(6);
  927. if (!all (all (p[ [1:3] ; ] == p[ [1:3]' ; ]))) { error (); }
  928. if (!all (all (p[ ; [1:3] ] == p[ ; [1:3]' ]))) { error (); }
  929. if (!all (all (p[ [2:6] ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
  930. if (!all (all (p[ [2:6]' ; [2:6] ] == p[ [2:6]'; [2:6]' ]))) { error (); }
  931. if (!all (all (p[ [6:12]' ] == p[ [6:12] ]))) { error (); }
  932.  
  933. //
  934. //  Sub-Matrix promotion
  935. //
  936.  
  937. if(m2[2;2] != 4) { error(); }
  938. if(any(m2[2;] != [3,4])) { error(); }
  939. if(any(m2[;2] != [2,4]')) { error(); }
  940. i=2;j=1;
  941. if(m2[i;j] != 3) { error(); }
  942. i=1;j=2;
  943. if(m2[i;j] != 2) { error(); }
  944. m = [1,2,3;4,5,6;7,8,9];
  945.  
  946. if(any(m[1;1,2] != [1,2])) 
  947. {
  948.   error();
  949. }
  950.  
  951. if(any(m[1,2;1] != [1;4])) 
  952. {
  953.   error();
  954. }
  955.  
  956. if(any (any (m[1,2;1,2] != [1,2;4,5]))) 
  957. {
  958.   error();
  959. }
  960.  
  961. //
  962.  
  963. if(m3[2;2] != (5+6i)) { error(); }
  964. if(any(m3[2;] != [3+4i,5+6i])) { error(); }
  965. if(any(m3[;2] != conj([2+3i,5+6i]'))) { error(); }
  966.  
  967. //
  968. //  Automatic creation, extension
  969. //
  970.  
  971. if(any (any ((mm[3;3]=10) != [0,0,0;0,0,0;0,0,10]))) { error(); }
  972. a=[1,2,3;4,5,6];
  973. if(any (any ((a[3;1]=10) != [1,2,3;4,5,6;10,0,0]))) { error(); }
  974. a=[1,2;3,4];
  975. if(any (any ((a[3,4;3,4]=[5,6;7,8]) != [1,2,0,0;3,4,0,0;0,0,5,6;0,0,7,8]))) 
  976. {
  977.   error();
  978. }
  979.  
  980. mmm[2;] = a[1;];
  981.  
  982. //
  983. //  Matrix binary operations
  984. //
  985.  
  986. a = m2; b = [5,6;7,8];
  987. if(any (any (a+a != [2,4;6,8]))) { error(); }
  988. if(any (any (a-a != zeros(2,2)))) { error(); }
  989. if(any (any (2+a != [3,4;5,6]))) { error(); }
  990. if(any (any (2-a != [1,0;-1,-2]))) { error(); }
  991. if(any (any (a-2 != [-1,0;1,2]))) { error(); }
  992. if(any (any (2*a != [2,4;6,8]))) { error(); }
  993. if(any (any ((a./2 != [0.5,1;1.5,2])))) { error(); }
  994. if(any (any (a*a != [7,10;15,22]))) { error(); }
  995. if(any (any (a*a*a != [37,54;81,118]))) { error(); }
  996. if(any (any (a .* a != [1,4;9,16]))) { error(); }
  997. if(any (any (a./a != [1,1;1,1]))) { error(); }
  998. if(any (any (a' != [1,3;2,4]))) { error(); }
  999.  
  1000. if(any(any(rand(3,3)^0 != eye(3,3)))) { error(); }
  1001. if(any(any(rand(3,3).^0 != ones(3,3)))) { error(); }
  1002. if(any(any(rand(1,3).^0 != ones(1,3)))) { error(); }
  1003. if(any(any(rand(3,1).^0 != ones(3,1)))) { error(); }
  1004. if(any(any(1.^zeros(3,1) != ones(3,1)))) { error(); }
  1005. if(any(any(1.^zeros(1,3) != ones(1,3)))) { error(); }
  1006.  
  1007. if(any ([1;2;3]+[4;5;6] != [5;7;9])) 
  1008. {
  1009.   error();
  1010. }
  1011.  
  1012. if(any ([1;2;3]-[4;5;6] != [-3;-3;-3])) 
  1013. {
  1014.   error();
  1015. }
  1016.  
  1017. if(any ([2;2;2] ./ [1;1;1] != [2;2;2])) 
  1018. {
  1019.   error();
  1020. }
  1021.  
  1022. if(any ([1;2;3] .* [4;5;6] != [4;10;18])) 
  1023. {
  1024.   error();
  1025. }
  1026.  
  1027. if (type (1^.33333) != "real") { error (); }
  1028. if (type (1^(1/3)) != "real") { error (); }
  1029. if (type ([1]^.33333) != "real") { error (); }
  1030. if (type (1^[.33333]) != "real") { error (); }
  1031. if (type ([1]^[.33333]) != "real") { error (); }
  1032.  
  1033. //
  1034. // Test row-wise matrix addition
  1035. //
  1036.  
  1037. a = [1,2,3]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1038. if (!all (all (a .+ b == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  1039. if (!all (all (b .+ a == [2,4,6;5,7,9;8,10,12;11,13,15]))) { error (); }
  1040.  
  1041. a = [1,2,3] + [1,2,3]*1i; 
  1042. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  1043. c = [2,4,6;5,7,9;8,10,12;11,13,15] + [2,4,6;5,7,9;8,10,12;11,13,15]*1i;
  1044. if (!all (all (a .+ b == c))) { error (); }
  1045. if (!all (all (b .+ a == c))) { error (); }
  1046.  
  1047. printf("\t\tpassed matrix row-wise add test...\n");
  1048.  
  1049. //
  1050. // Test row-wise matrix subtraction
  1051. //
  1052.  
  1053. a = [1,1,1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1054. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1055. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1056.  
  1057. a = [1,1,1] + [1,1,1]*1i;
  1058. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  1059. c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
  1060. if (!all (all (a .- b == -c))) { error (); }
  1061. if (!all (all (b .- a ==  c))) { error (); }
  1062.  
  1063. printf("\t\tpassed matrix row-wise subtraction test...\n");
  1064.  
  1065. //
  1066. // Test col-wise matrix addition
  1067. //
  1068.  
  1069. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1070. if (!all (all (a .+ b == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  1071. if (!all (all (b .+ a == [2,3,4;5,6,7;8,9,10;11,12,13]))) { error (); }
  1072.  
  1073. a = [1;1;1;1] + [1;1;1;1]*1i;
  1074. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  1075. c = [2,3,4;5,6,7;8,9,10;11,12,13] + [2,3,4;5,6,7;8,9,10;11,12,13]*1i;
  1076. if (!all (all (a .+ b == c))) { error (); }
  1077. if (!all (all (b .+ a == c))) { error (); }
  1078.  
  1079. printf("\t\tpassed matrix col-wise add test...\n");
  1080.  
  1081. //
  1082. // Test col-wise matrix subtraction
  1083. //
  1084.  
  1085. a = [1;1;1;1]; b = [1,2,3;4,5,6;7,8,9;10,11,12];
  1086. if (!all (all (a .- b == -[0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1087. if (!all (all (b .- a ==  [0,1,2;3,4,5;6,7,8;9,10,11]))) { error (); }
  1088.  
  1089. a = [1;1;1;1] + [1;1;1;1]*1i;
  1090. b = [1,2,3;4,5,6;7,8,9;10,11,12] + [1,2,3;4,5,6;7,8,9;10,11,12]*1i;
  1091. c = [0,1,2;3,4,5;6,7,8;9,10,11] + [0,1,2;3,4,5;6,7,8;9,10,11]*1i;
  1092. if (!all (all (a .- b == -c))) { error (); }
  1093. if (!all (all (b .- a ==  c))) { error (); }
  1094.  
  1095. printf("\t\tpassed matrix col-wise subtraction test...\n");
  1096.  
  1097. a = [1,2,3];
  1098. b = [1,2,3;4,5,6;7,8,9];
  1099. c = [1,4,9;4,10,18;7,16,27];
  1100.  
  1101. if (!all (all (a .* b == c))) { error (); }
  1102. if (!all (all (b .* a == c))) { error (); }
  1103.  
  1104. za = a + rand (size (a))*1j;
  1105. zb = b + rand (size (b))*1j;
  1106.  
  1107. if (!all (all (za .* zb == [za;za;za] .* zb))) { error (); }
  1108. if (!all (all (zb .* za == zb .* [za;za;za]))) { error (); }
  1109. if (!all (all (a .* zb == [a;a;a] .* zb))) { error (); }
  1110. if (!all (all (zb .* a == zb .* [a;a;a]))) { error (); }
  1111. if (!all (all (za .* b == [za;za;za] .* b))) { error (); }
  1112. if (!all (all (b .* za == b .* [za;za;za]))) { error (); }
  1113.  
  1114. printf("\t\tpassed matrix row-wise multiplication test...\n");
  1115.  
  1116. a = [1,2,3];
  1117. b = [1,2,3;4,6,6;7,8,9];
  1118. c = [1,1,1;4,3,2;7,4,3];
  1119.  
  1120. if (!all (all (b ./ a == c))) { error (); }
  1121. if (!all (all ([a;a;a] ./ b == a ./ b))) { error (); }
  1122. if (!all (all (b ./ [a;a;a] == b ./ a))) { error (); }
  1123.  
  1124. za = a + rand (size (a))*1j;
  1125. zb = b + rand (size (b))*1j;
  1126.  
  1127. if (!all (all ([za;za;za] ./ zb == za ./ zb))) { error (); }
  1128. if (!all (all (zb ./ [za;za;za] == zb ./ za))) { error (); }
  1129. if (!all (all ([a;a;a] ./ zb == a ./ zb))) { error (); }
  1130. if (!all (all (zb ./ [a;a;a] == zb ./ a))) { error (); }
  1131. if (!all (all ([za;za;za] ./ b == za ./ b))) { error (); }
  1132. if (!all (all (b ./ [za;za;za] == b ./ za))) { error (); }
  1133.  
  1134. printf("\t\tpassed matrix row-wise division test...\n");
  1135.  
  1136. a = [1;2;3];
  1137. b = [1,2,3;4,5,6;7,8,9];
  1138.  
  1139. if (!all (all (a .* b == [a,a,a] .* b))) { error (); }
  1140. if (!all (all (b .* a == b .* [a,a,a]))) { error (); }
  1141.  
  1142. za = a + rand (size (a))*1j;
  1143. zb = b + rand (size (b))*1j;
  1144.  
  1145. if (!all (all (za .* zb == [za,za,za] .* zb))) { error (); }
  1146. if (!all (all (zb .* za == zb .* [za,za,za]))) { error (); }
  1147. if (!all (all (za .* b == [za,za,za] .* b))) { error (); }
  1148. if (!all (all (b .* za == b .* [za,za,za]))) { error (); }
  1149. if (!all (all (a .* zb == [a,a,a] .* zb))) { error (); }
  1150. if (!all (all (zb .* a == zb .* [a,a,a]))) { error (); }
  1151.  
  1152. printf("\t\tpassed matrix column-wise multiplication test...\n");
  1153.  
  1154. a = [1;2;3];
  1155. b = [1,2,3;4,6,6;7,8,9];
  1156.  
  1157. if (!all (all ([a,a,a] ./ b == a ./ b))) { error (); }
  1158. if (!all (all (b ./ [a,a,a] == b ./ a))) { error (); }
  1159.  
  1160. za = a + rand (size(a))*1j;
  1161. zb = b + rand (size(b))*1j;
  1162.  
  1163. if (!all (all ([za,za,za] ./ zb == za ./ zb))) { error (); }
  1164. if (!all (all (zb ./ [za,za,za] == zb ./ za))) { error (); }
  1165. if (!all (all ([za,za,za] ./ b == za ./ b))) { error (); }
  1166. if (!all (all (b ./ [za,za,za] == b ./ za))) { error (); }
  1167. if (!all (all ([a,a,a] ./ zb == a ./ zb))) { error (); }
  1168. if (!all (all (zb ./ [a,a,a] == zb ./ a))) { error (); }
  1169.  
  1170. printf("\t\tpassed matrix column-wise division test...\n");
  1171.  
  1172.  
  1173. //
  1174. //--------------------- Test COMPLEX MATRICES --------------------------
  1175. //
  1176. //  Automatic creation, extension
  1177. //
  1178. printf("\t\tcomplex-matrices\n");
  1179. if(any (any ((mm[3;3]=10+10i) != [0,0,0;0,0,0;0,0,10+10i]))) { error(); }
  1180. a=[1,2,3;4,5,6];
  1181. if(any (any ((a[3;1]=10+10i) != [1,2,3;4,5,6;10+10i,0,0]))) { error(); }
  1182. //
  1183. a = m3;
  1184. if(any (any (a+a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  1185. if(any (any (a-a != zeros(2,2)))) { error(); }
  1186. if(any (any (2+a != [3+2i,4+3i;5+4i,7+6i]))) { error(); }
  1187. if(any (any (2-a != [1-2i,0-3i;-1-4i,-3-6i]))) { error(); }
  1188. if(any (any (a-2 != [-1+2i,0+3i;1+4i,3+6i]))) { error(); }
  1189. if(any (any (2*a != [2+4i,4+6i;6+8i,10+12i]))) { error(); }
  1190. if(any (any (a./2 != [.5+1i,1+1.5i;1.5+2i,2.5+3i]))) { error(); }
  1191. if(any (any (a*a != [-9+21i,-12+34i;-14+48i,-17+77i]))) { error(); }
  1192. if(any (any (a*a*a != [-223+57i,-345+113i;-469+183i,-719+337i]))) { error(); }
  1193. if(any (any (a .* a != [-3+4i,-5+12i;-7+24i,-11+60i]))) { error(); }
  1194. //
  1195. // The following test may not work on some machines
  1196. //
  1197. if(any (any (a./a != [1,1;1,1]))) { 
  1198.   printf("\t\t***complex division inaccuracy, check manually***\n");
  1199. }
  1200.  
  1201. if(any (any (a' != conj([1+2i,3+4i;2+3i,5+6i])))) { error(); }
  1202. //  
  1203. //--------------------- Test NULL MATRICES -------------------------
  1204. //
  1205. printf("\t\tnull-matrices\n");
  1206. // Create a NULL matrix
  1207. mnull = [];
  1208. // Test it with SCALARS
  1209. if( any([1,mnull] != 1)) {
  1210.   error();
  1211. }
  1212. if( any([mnull,1] != 1)) {
  1213.   error();
  1214. }
  1215. // Test with MATRIX construction
  1216. m = [1,2;3,4;5,6];
  1217. if( any([mnull;1] != [1])) {
  1218.   error();
  1219. }
  1220. if( any([1;mnull] != [1])) {
  1221.   error();
  1222. }
  1223. if( any([mnull;1,2,3] != [1,2,3])) {
  1224.   error();
  1225. }
  1226. if( any([1,2,3;mnull] != [1,2,3])) {
  1227.   error();
  1228. }
  1229. if(any (any ([mnull,m] != m))) {
  1230.   error();
  1231. }
  1232. if(any (any ([m,mnull] != m))) {
  1233.   error();
  1234. }
  1235. if(any (any ([mnull;m] != m))) {
  1236.   error();
  1237. }
  1238. if(any (any ([m;mnull] != m))) {
  1239.   error();
  1240.  
  1241. mnull = matrix();
  1242. mnull[1] = [1];
  1243. }
  1244.  
  1245. //--------------------- Test Matrix Multiply --------------------------
  1246.  
  1247. i = sqrt(-1);
  1248. a = [1,2,3;4,5,6;7,8,9];
  1249. b = [4,5,6;7,8,9;10,11,12];
  1250. c = [ 48,  54,  60 ;
  1251.      111, 126, 141 ;
  1252.      174, 198, 222 ];
  1253.  
  1254. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  1255.  
  1256. az = a + b*i;
  1257. bz = b + a*i;
  1258.  
  1259. cz = [-18+141*i , -27+162*i , -36+183*i ;
  1260.         9+240*i ,   0+279*i ,  -9+318*i ;
  1261.        36+339*i ,  27+396*i ,  18+453*i ];
  1262.  
  1263. czz = [ 48+30*i ,  54+36*i  ,  60+42*i ;
  1264.        111+66*i , 126+81*i  , 141+96*i ;
  1265.        174+102*i, 198+126*i , 222+150*i ];
  1266.  
  1267. czzz = [ 48+111*i ,  54+126*i ,  60+141*i ;
  1268.         111+174*i , 126+198*i , 141+222*i ;
  1269.         174+237*i , 198+270*i , 222+303*i ];
  1270.  
  1271. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  1272. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  1273. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  1274.  
  1275. a = [a,a];
  1276. b = [b;b];
  1277. c = [  96 , 108 , 120 ;
  1278.       222 , 252 , 282 ;
  1279.       348 , 396 , 444 ];
  1280.  
  1281. if (any (any (c != a*b))) { error ("failed Real-Real Multiply"); }
  1282.  
  1283. az = [az,az];
  1284. bz = [bz;bz];
  1285.  
  1286. cz = [  -36+282*i ,  -54+324*i ,  -72+366*i ;
  1287.          18+480*i ,    0+558*i ,  -18+636*i ;
  1288.          72+678*i ,   54+792*i ,   36+906*i ];
  1289.  
  1290. czz = [  96+60*i  , 108+72*i  , 120+84*i  ;
  1291.         222+132*i , 252+162*i , 282+192*i ;
  1292.         348+204*i , 396+252*i , 444+300*i ];
  1293.  
  1294. czzz = [  96+222*i , 108+252*i , 120+282*i ;
  1295.          222+348*i , 252+396*i , 282+444*i ;
  1296.          348+474*i , 396+540*i , 444+606*i ];
  1297.  
  1298. if (any (any (cz != az*bz)))  { error ("failed Complex-Complex Multiply"); }
  1299. if (any (any (czz != a*bz)))  { error ("failed Real-Complex Multiply"); }
  1300. if (any (any (czzz != az*b))) { error ("failed Complex-Real Multiply"); }
  1301.  
  1302. printf("\t\tpassed matrix multiply test...\n");
  1303.  
  1304. //--------------------- Test STRING MATRICES --------------------------
  1305. //
  1306. printf("\t\tstring-matrices\n");
  1307. sm = ["s1","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"];
  1308. if(sm[1] != "s1") { error(); }
  1309. if( sm[1;3] != "sm3" ) { error(); }
  1310. if(any(sm[2,3;3] != ["xxx3";"yyy3"]) ) { error(); }
  1311. if(any (any ((sm[1;1]="xx")!=["xx","sm2","sm3";"x1","x2","xxx3";"y1","y2","yyy3"])))
  1312. {
  1313.   error();
  1314. }
  1315. if( ((strm[1] = "strm")[1]) != "strm" ) { error(); }
  1316.  
  1317. //  Test string-matrix add functionality
  1318.  
  1319. sm = sm[1,2;1,2];
  1320. if (any (any (("1_"+sm+"_2") != ["1_xx_2","1_sm2_2";"1_x1_2","1_x2_2"]))) {error();}
  1321.  
  1322. if ("c"+["1"] != "c1") { error (); }
  1323. if (["c"]+"1" != "c1") { error (); }
  1324.  
  1325. printf("\tpassed matrix test...\n");
  1326.  
  1327. //
  1328. //---------------------------- List Tests --------------------------
  1329. //
  1330. //  List creation
  1331. listest = << << 11; 12 >>; << 21; 22>> >>;
  1332. if( listest.[1].[2] != 12 ) { error(); }
  1333. if(any(<<a=10;b=1:4;c=[1,2,3;4,5,6;7,8,9]>>.b != [1,2,3,4])) { error(); }
  1334. mlist.[0] = m;
  1335. if(any(any(mlist.[0] != m))) { error(); }
  1336. printf("\tpassed list test...\n");
  1337.  
  1338. // Test list functions...
  1339.  
  1340. listtest.fun = function ( a ) { return sin(a); };
  1341. if (!(listtest.fun(0.5) == sin(0.5))) { error() ; }
  1342.  
  1343. //
  1344. // Reset random number generator seed...
  1345. //
  1346.  
  1347. rand("default");
  1348.  
  1349. //
  1350. //-------------------------- Test printf () --------------------------
  1351. //
  1352.  
  1353. sprintf (tmp, "%*.*d %*.*d %s %*.*f f\n", 5,3,2, 8,7,3, "string", 3, 4, 1234e-2);
  1354. if (!(tmp == "  002  0000003 string 12.3400 f\n")) { error ("sprintf() error"); }
  1355. printf("\tpassed sprintf test...\n");
  1356.  
  1357. //
  1358. //------------------------- Test strtod()  ----------------------------
  1359. //
  1360.  
  1361. if (123.456 != strtod ("123.456")) { error (); }
  1362. if (!all (all ([1,2;3,4] == strtod (["1","2";"3","4"]))))
  1363.   { error (); }
  1364. printf("\tpassed strtod test...\n");
  1365.  
  1366. //
  1367. //------------------------- Test getline()  ---------------------------
  1368. //
  1369. //
  1370.  
  1371. close( "test.getline" );
  1372.  
  1373. x = getline( "test.getline" );
  1374. if ( x.[1] !=  123.456 ) { error(); }
  1375. if ( x.[2] != -123.456 ) { error(); }
  1376. if ( x.[3] !=  123.456 ) { error(); }
  1377.  
  1378. x = getline( "test.getline" );
  1379. if ( x.[1] !=  .123 ) { error(); }
  1380. if ( x.[2] != -.123 ) { error(); }
  1381. if ( x.[3] !=  .123 ) { error(); }
  1382.  
  1383. x = getline( "test.getline" );
  1384. if ( x.[1] !=  123 ) { error(); }
  1385. if ( x.[2] != -123 ) { error(); }
  1386. if ( x.[3] !=  123 ) { error(); }
  1387.  
  1388. x = getline( "test.getline" );
  1389. if ( x.[1] !=  1e6 ) { error(); }
  1390. if ( x.[2] != -1e6 ) { error(); }
  1391. if ( x.[3] !=  1e6 ) { error(); }
  1392.  
  1393. x = getline( "test.getline" );
  1394. if ( x.[1] !=  1e5 ) { error(); }
  1395. if ( x.[2] != -1e5 ) { error(); }
  1396. if ( x.[3] !=  1e5 ) { error(); }
  1397.  
  1398. x = getline( "test.getline" );
  1399. if ( x.[1] !=  123.456e3 ) { error(); }
  1400. if ( x.[2] != -123.456e3 ) { error(); }
  1401. if ( x.[3] !=  123.456e3 ) { error(); }
  1402.  
  1403. x = getline( "test.getline" );
  1404. if ( x.[1] !=  123.456e3 ) { error(); }
  1405. if ( x.[2] != -123.456e3 ) { error(); }
  1406. if ( x.[3] !=  123.456e3 ) { error(); }
  1407.  
  1408. x = getline( "test.getline" );
  1409. if ( x.[1] !=  123.456e-3 ) { error(); }
  1410. if ( x.[2] != -123.456e-3 ) { error(); }
  1411. if ( x.[3] !=  123.456e-3 ) { error(); }
  1412.  
  1413. x = getline( "test.getline" );
  1414. if ( x.[1] !=  .123e3 ) { error(); }
  1415. if ( x.[2] != -.123e3 ) { error(); }
  1416. if ( x.[3] !=  .123e3 ) { error(); }
  1417.  
  1418. x = getline( "test.getline" );
  1419. if ( x.[1] !=  123e3 ) { error(); }
  1420. if ( x.[2] != -123e3 ) { error(); }
  1421. if ( x.[3] !=  123e3 ) { error(); }
  1422.  
  1423. x = getline( "test.getline" );
  1424. if ( x.[1] != "123abc" ) { error(); }
  1425. if ( x.[2] != "abc123" ) { error(); }
  1426. if ( x.[3] != "123.abc" ) { error(); }
  1427.  
  1428. x = getline( "test.getline" );
  1429. if ( x.[1] != "quoted string" ) { error(); }
  1430. if ( x.[2] != "q string with escapes \n \t \" " ) { error(); }
  1431.  
  1432. x = getline( "test.getline" );
  1433. if ( x.[1] != "quoted string" ) { error(); }
  1434. if ( x.[2] !=  1.23e3 ) { error(); }
  1435. if ( x.[3] !=  100 ) { error(); }
  1436. if ( x.[4] != "q string with escapes \n \t \" " ) { error(); }
  1437. if ( x.[5] !=  200 ) { error(); }
  1438.  
  1439. printf("\tpassed getline() test...\n");
  1440.  
  1441. //
  1442. //---------------------- Test readb()/writeb()  --------------------
  1443. //
  1444.  
  1445. a = rand (5,5);
  1446. z = rand(3,3) + rand(3,3)*1j;
  1447. strm = what()[1:5;1:5];
  1448. pi = 4*atan(1.0);
  1449. sc = 2*pi;
  1450. str = "this is a sample string\ttab";
  1451. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1452.  
  1453. writeb ("jnk.rb", a, z, strm, sc, str, l);
  1454.  
  1455. # Set aside matrices for later tests
  1456.  
  1457. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1458.  
  1459. clear (a, z, strm, sc, str, l);
  1460.  
  1461. close ("jnk.rb");
  1462.  
  1463. readb ("jnk.rb");
  1464.  
  1465. #
  1466. # Now do checks
  1467. #
  1468.  
  1469. if (!all (all (a == check.a))) { error (); }
  1470. if (!all (all (z == check.z))) { error (); }
  1471. if (!all (all (strm == check.strm))) { error (); }
  1472. if (sc != check.sc) { error (); }
  1473. if (str != check.str) { error (); }
  1474.  
  1475. if (length (l) != 5) { error (); }
  1476. if (!all (all (l.a == check.a))) { error (); }
  1477. if (!all (all (l.z == check.z))) { error (); }
  1478. if (!all (all (l.strm == check.strm))) { error (); }
  1479. if (l.sc != check.sc) { error (); }
  1480. if (l.str != check.str) { error (); }
  1481.  
  1482.  
  1483. printf("\tpassed binary I/O test...\n");
  1484.  
  1485. //
  1486. //---------------------- Test read()/write()  --------------------
  1487. //
  1488.  
  1489. a = rand (5,5);
  1490. z = rand(3,3) + rand(3,3)*1j;
  1491. strm = what()[1:5;1:5];
  1492. pi = 4*atan(1.0);
  1493. sc = 2*pi;
  1494. str = "this is a sample string\ttab";
  1495. l = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1496.  
  1497. write ("jnk.ra", a, z, strm, sc, str, l);
  1498.  
  1499. # Set aside matrices for later tests
  1500.  
  1501. check = <<a = a; z = z; strm = strm; sc = sc; str = str>>;
  1502.  
  1503. clear (a, z, strm, sc, str, l);
  1504.  
  1505. close ("jnk.ra");
  1506.  
  1507. read ("jnk.ra");
  1508.  
  1509. #
  1510. # Now do checks
  1511. #
  1512.  
  1513. if (a.nr != 5 || a.nc != 5) { error (); }
  1514. if (z.nr != 3 || z.nc != 3) { error (); }
  1515. if (strm.nr != 5 || strm.nc != 5) { error (); }
  1516. if (str != check.str) { error (); }
  1517.  
  1518. if (length (l) != 5) { error (); }
  1519. if (l.a.nr != 5 || l.a.nc != 5) { error (); }
  1520. if (l.z.nr != 3 || l.z.nc != 3) { error (); }
  1521. if (l.strm.nr != 5 || l.strm.nc != 5) { error (); }
  1522. if (l.str != check.str) { error (); }
  1523.  
  1524. printf("\tpassed ascii I/O test...\n");
  1525.  
  1526. //
  1527. //-------------------------- Test eval () ------------------------------
  1528. //
  1529.  
  1530. if (1 + 2 != eval("1 + 2")) { error ("eval() error"); }
  1531. x = function (s, a) { return eval(s); };
  1532. str = "yy = 2 + x(\"2*a\",3)";
  1533. z = eval(str);
  1534. if (z != yy) { error ("eval() error"); }
  1535. printf("\tpassed eval test...\n");
  1536.  
  1537.  
  1538. //-------------------------- ------------------------------
  1539.  
  1540. printf ("\ttest builtin function correctness...\n");
  1541.  
  1542. //-------------------------- ------------------------------
  1543.  
  1544. //
  1545. //-------------------------- Test abs () ------------------------------
  1546. //
  1547.  
  1548. A = rand(5,5);
  1549. T = ( A == abs (A));
  1550. if (!all (all (A))) { error ("abs() incorrect"); }
  1551. printf("\tabs()");
  1552.  
  1553.  
  1554. //
  1555. //-------------------------- Test max () ------------------------------
  1556. //
  1557.  
  1558. A = [1,10,100;2,20,200;1,2,3];
  1559. B = A./2;
  1560. ZA = A + rand (3,3)*A*1i;
  1561. ZB = B + rand (3,3)*B*1i;
  1562. if (!all (max (A) == [2,20,200])) { error( "max() incorrect"); }
  1563. if (max (max(A)) != 200) { error ("max() incorrect"); }
  1564. if (any (any (max (A, B) != A))) { error (); }
  1565. if (any (any (max (B, A) != max (A, B)))) { error (); }
  1566. if (any (any (max (ZB, ZA) != max (ZA, ZB)))) { error (); }
  1567. if (any (any (max (B, ZA) != max (ZA, B)))) { error (); }
  1568. if (any (any (max (ZB, A) != max (A, ZB)))) { error (); }
  1569. printf("\tmax()");
  1570.  
  1571. //
  1572. //-------------------------- Test min () ------------------------------
  1573. //
  1574.  
  1575. if (!all (min (A) == [1,2,3])) { error( "min() incorrect"); }
  1576. if (min (min(A)) != 1) { error ("min() incorrect"); }
  1577. if (any (any (min (A, B) != B))) { error (); }
  1578. if (any (any (min (B, A) != min (A, B)))) { error (); }
  1579. if (any (any (min (ZB, ZA) != min (ZA, ZB)))) { error (); }
  1580. if (any (any (min (B, ZA) != min (ZA, B)))) { error (); }
  1581. if (any (any (min (ZB, A) != min (A, ZB)))) { error (); }
  1582. printf("\tmin()");
  1583.  
  1584. //
  1585. //-------------------------- Test maxi () -----------------------------
  1586. //
  1587.  
  1588. if (!all (maxi (A) == [2,2,2])) { error( "maxi() incorrect"); }
  1589. if (maxi (maxi(A)) != 1) { error ("maxi() incorrect"); }
  1590. printf("\tmaxi()");
  1591.  
  1592. //
  1593. //-------------------------- Test mini () -----------------------------
  1594. //
  1595.  
  1596. if (!all (mini (A) == [1,3,3])) { error( "mini() incorrect"); }
  1597. if (mini (mini(A)) != 1) { error ("mini() incorrect"); }
  1598. printf("\tmini()");
  1599.  
  1600. //
  1601. //-------------------------- Test floor () ----------------------------
  1602. //
  1603.  
  1604. if (floor (1.9999) != 1) { error ("floor() output incorrect"); }
  1605. if (!all (all (floor ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  1606.   error ("floor output incorrect");
  1607. }
  1608. printf("\tfloor()");
  1609.  
  1610. //
  1611. //-------------------------- Test ceil () 0----------------------------
  1612. //
  1613.  
  1614. if (ceil (1.9999) != 2) { error ("ceil() output incorrect"); }
  1615. if (!all (all (ceil ([1.99,1.99;2.99,2.99]) == [2,2;3,3]))) {
  1616.   error ("ceil output incorrect");
  1617. }
  1618. printf("\tceil()\n");
  1619.  
  1620. //
  1621. //-------------------------- Test round () ----------------------------
  1622. //
  1623.  
  1624. if (round (1.8) != 2) { error ("round() output incorrect"); }
  1625. if (round (1.4) != 1) { error ("round() output incorrect"); }
  1626. if (!all (all (round ([1.99,1.99;2.4,2.4]) == [2,2;2,2]))) {
  1627.   error ("round output incorrect");
  1628. }
  1629. printf("\tround()");
  1630.  
  1631. //
  1632. //-------------------------- Test sum () ------------------------------
  1633. //
  1634.  
  1635. S = [1:4; 4:7; 8:11];
  1636. if (sum (S[1;]) != 10) { error ("sum() incorrect"); }
  1637. if (sum (S[3;]) != 38) { error ("sum() incorrect"); }
  1638. if (!all (all (sum (S) == [13,16,19,22]))) { error ("sum() incorrect"); }
  1639. printf("\tsum()");
  1640.  
  1641. //
  1642. //-------------------------- Test int () ------------------------------
  1643. //
  1644.  
  1645. if (int (1.9999) != 1) { error ("int() output incorrect"); }
  1646. if (!all (all (int ([1.99,1.99;2.99,2.99]) == [1,1;2,2]))) {
  1647.   error ("int() output incorrect");
  1648. }
  1649. printf("\tint()");
  1650.  
  1651. //
  1652. //-------------------------- Test mod () ------------------------------
  1653. //
  1654.  
  1655. if (mod (1,1) != 0) { error ("mod() output incorrect"); }
  1656. if (mod (4,2) != 0) { error ("mod() output incorrect"); }
  1657. if (mod (3,2) != 1) { error ("mod() output incorrect"); }
  1658. if (mod (5,3) != 2) { error ("mod() output incorrect"); }
  1659. printf("\tmod()");
  1660.  
  1661. //
  1662. //-------------------------- Test find () ------------------------------
  1663. //
  1664.  
  1665. if (find ([0,1]) != 2) { error ("find() output incorrect"); }
  1666. if (find ([1,0]) != 1) { error ("find() output incorrect"); }
  1667. if (find ([0,1+1i]) != 2) { error ("find() output incorrect"); }
  1668. if (find ([1+0i,0]) != 1) { error ("find() output incorrect"); }
  1669. if (find ([0+1i,0]) != 1) { error ("find() output incorrect"); }
  1670.  
  1671. printf("\tfind()\n");
  1672.  
  1673.  
  1674. //-------------------------- ------------------------------
  1675.  
  1676. printf ("\ttest builtin function operation...\n");
  1677.  
  1678. //-------------------------- ------------------------------
  1679.  
  1680. //
  1681. // At least use most of the builtins....
  1682. //
  1683.  
  1684. A = rand(5,5);
  1685. Z = rand(5,5) + rand(5,5)*1j;
  1686. S = what()[3:6;2:4];
  1687.  
  1688. abs (A);
  1689. abs ([A]);
  1690. abs (3.14);
  1691.  
  1692. abs (Z);
  1693. abs ([Z]);
  1694. abs (3.14j);
  1695.  
  1696. printf ("\tabs");
  1697.  
  1698. acos (A);
  1699. acos ([A]);
  1700. acos (3.14);
  1701.  
  1702. acos (Z);
  1703. acos ([Z]);
  1704. acos (3.14j);
  1705.  
  1706. printf ("\tacos");
  1707.  
  1708. all (A);
  1709. all ([A]);
  1710. all (3.14);
  1711.  
  1712. all (Z);
  1713. all ([Z]);
  1714. all (3.14j);
  1715.  
  1716. printf ("\tall");
  1717.  
  1718. any (A);
  1719. any ([A]);
  1720. any (3.14);
  1721.  
  1722. any (Z);
  1723. any ([Z]);
  1724. any (3.14j);
  1725.  
  1726. printf ("\tany");
  1727.  
  1728. asin (A);
  1729. asin ([A]);
  1730. asin (3.14);
  1731.  
  1732. asin (Z);
  1733. asin ([Z]);
  1734. asin (3.14j);
  1735.  
  1736. printf ("\tasin");
  1737.  
  1738. atan (A);
  1739. atan ([A]);
  1740. atan (3.14);
  1741.  
  1742. atan (Z);
  1743. atan ([Z]);
  1744. atan (3.14j);
  1745.  
  1746. printf ("\tatan");
  1747.  
  1748. atan2 (A,A);
  1749. atan2 ([A],[A]);
  1750. atan2 (3.14,3.14);
  1751.  
  1752. #atan2 (Z,Z);
  1753. #atan2 ([Z],[Z]);
  1754. #atan2 (3.14j,3.14j);
  1755.  
  1756. printf ("\tatan2");
  1757. printf ("\n");
  1758.  
  1759. cd (":");
  1760. printf ("\tcd");
  1761.  
  1762. ceil (A);
  1763. ceil ([A]);
  1764. ceil (3.14);
  1765.  
  1766. ceil (Z);
  1767. ceil ([Z]);
  1768. ceil (3.14j);
  1769.  
  1770. printf ("\tceil");
  1771.  
  1772. class (A);
  1773. class ([A]);
  1774. class (3.14);
  1775. class ("string");
  1776. class (S);
  1777. class ([S]);
  1778. class (<< rand(3,3); rand(4,4) >>);
  1779.  
  1780. class (Z);
  1781. class ([Z]);
  1782. class (3.14j);
  1783.  
  1784. printf ("\tclass");
  1785.  
  1786. clear (rand(3,3));
  1787. clear ([rand(3,3)]);
  1788.  
  1789. printf ("\tclear");
  1790.  
  1791. conj (A);
  1792. conj ([A]);
  1793. conj (3.14);
  1794.  
  1795. conj (Z);
  1796. conj ([Z]);
  1797. conj (3.14j);
  1798.  
  1799. printf ("\tconj");
  1800.  
  1801. cos (A);
  1802. cos ([A]);
  1803. cos (3.14);
  1804.  
  1805. cos (Z);
  1806. cos ([Z]);
  1807. cos (3.14j);
  1808.  
  1809. printf ("\tcos");
  1810.  
  1811. diag (A);
  1812. diag ([A]);
  1813. diag (3.14);
  1814.  
  1815. diag (Z);
  1816. diag ([Z]);
  1817. diag (3.14j);
  1818.  
  1819. printf ("\tdiag");
  1820. printf ("\n");
  1821.  
  1822. exist (A);
  1823. exist (Z);
  1824.  
  1825. printf ("\texist");
  1826.  
  1827. exp (A);
  1828. exp ([A]);
  1829. exp (3.14);
  1830.  
  1831. exp (Z);
  1832. exp ([Z]);
  1833. exp (3.14j);
  1834.  
  1835. printf ("\texp");
  1836.  
  1837. find (A);
  1838. find ([A]);
  1839. find (3.14);
  1840.  
  1841. find (Z);
  1842. find ([Z]);
  1843. find (3.14j);
  1844.  
  1845. printf ("\tfind");
  1846.  
  1847. floor (A);
  1848. floor ([A]);
  1849. floor (3.14);
  1850.  
  1851. floor (Z);
  1852. floor ([Z]);
  1853. floor (3.14j);
  1854.  
  1855. printf ("\tfloor");
  1856.  
  1857. format (1,1);
  1858. format ([2],[3]);
  1859. format ();
  1860. format ([3], 3);
  1861. format ();
  1862.  
  1863. printf ("\tformat");
  1864.  
  1865. imag (A);
  1866. imag ([A]);
  1867. imag (3.14);
  1868.  
  1869. imag (Z);
  1870. imag ([Z]);
  1871. imag (3.14j);
  1872.  
  1873. printf ("\timag");
  1874.  
  1875. inf ();
  1876. printf ("\tinf");
  1877. printf ("\n");
  1878.  
  1879. int (A);
  1880. int ([A]);
  1881. int (3.14);
  1882.  
  1883. int (Z);
  1884. int ([Z]);
  1885. int (3.14j);
  1886.  
  1887. printf ("\tint");
  1888.  
  1889. issymm (A);
  1890. issymm ([A]);
  1891. issymm (3.14);
  1892.  
  1893. issymm (Z);
  1894. issymm ([Z]);
  1895. issymm (3.14j);
  1896.  
  1897. printf ("\tissymm");
  1898.  
  1899. length (A);
  1900. length ([A]);
  1901. length (3.14);
  1902.  
  1903. length (Z);
  1904. length ([Z]);
  1905. length (3.14j);
  1906.  
  1907. length (3.14);
  1908. length (S);
  1909. length (<< rand(2,2); rand(10,10); 3.14 >>);
  1910.  
  1911. printf ("\tlength");
  1912.  
  1913. log (A);
  1914. log ([A]);
  1915. log (3.14);
  1916.  
  1917. log (Z);
  1918. log ([Z]);
  1919. log (3.14j);
  1920.  
  1921. printf ("\tlog");
  1922.  
  1923. log10 (A);
  1924. log10 ([A]);
  1925. log10 (3.14);
  1926.  
  1927. log10 (Z);
  1928. log10 ([Z]);
  1929. log10 (3.14j);
  1930.  
  1931. printf ("\tlog10");
  1932.  
  1933. max (A);
  1934. max ([A]);
  1935. max (3.14);
  1936.  
  1937. max (Z);
  1938. max ([Z]);
  1939. max (3.14j);
  1940.  
  1941. printf ("\tmax");
  1942.  
  1943. maxi (A);
  1944. maxi ([A]);
  1945. maxi (3.14);
  1946.  
  1947. maxi (Z);
  1948. maxi ([Z]);
  1949. maxi (3.14j);
  1950.  
  1951. printf ("\tmaxi");
  1952. printf ("\n");
  1953.  
  1954. members (<< rand(2,2); rand(); rand(3,3) >>);
  1955.  
  1956. printf ("\tmembers");
  1957.  
  1958. min (A);
  1959. min ([A]);
  1960. min (3.14);
  1961.  
  1962. min (Z);
  1963. min ([Z]);
  1964. min (3.14j);
  1965.  
  1966. printf ("\tmin");
  1967.  
  1968. mini (A);
  1969. mini ([A]);
  1970. mini (3.14);
  1971.  
  1972. mini (Z);
  1973. mini ([Z]);
  1974. mini (3.14j);
  1975.  
  1976. printf ("\tmini");
  1977.  
  1978. mod (A,A);
  1979. mod ([A],A);
  1980. mod (3.14,2);
  1981.  
  1982. mod (Z,Z);
  1983. mod ([Z],Z);
  1984. mod (3.14j,2);
  1985.  
  1986. printf ("\tmod");
  1987.  
  1988. nan ();
  1989.  
  1990. ones(3,3);
  1991. ones([4,4]);
  1992.  
  1993. printf ("\tones");
  1994.  
  1995. prod (A);
  1996. prod ([A]);
  1997. prod (3.14);
  1998.  
  1999. prod (Z);
  2000. prod ([Z]);
  2001. prod (3.14j);
  2002.  
  2003. printf ("\tprod");
  2004.  
  2005. real (A);
  2006. real ([A]);
  2007. real (3.14);
  2008.  
  2009. real (Z);
  2010. real ([Z]);
  2011. real (3.14j);
  2012.  
  2013. printf ("\treal");
  2014. printf ("\n");
  2015.  
  2016. reshape (A, A.nr*A.nc, 1);
  2017. reshape (A, [A.nr*A.nc], [1]);
  2018. reshape (Z, Z.nr*Z.nc, 1);
  2019. reshape (Z, [Z.nr*Z.nc], [1]);
  2020.  
  2021. printf ("\treshape");
  2022.  
  2023. round (A);
  2024. round ([A]);
  2025. round (3.14);
  2026.  
  2027. round (Z);
  2028. round ([Z]);
  2029. round (3.14j);
  2030.  
  2031. printf ("\tround");
  2032.  
  2033. show(A);
  2034. show(Z);
  2035. show([A]);
  2036. show([Z]);
  2037.  
  2038. printf ("\tshow");
  2039.  
  2040. sign (A);
  2041. sign ([A]);
  2042. sign (3.14);
  2043.  
  2044. sign (Z);
  2045. sign ([Z]);
  2046. sign (3.14j);
  2047.  
  2048. printf ("\tsign");
  2049.  
  2050. sin (A);
  2051. sin ([A]);
  2052. sin (3.14);
  2053.  
  2054. sin (Z);
  2055. sin ([Z]);
  2056. sin (3.14j);
  2057.  
  2058. printf ("\tsin");
  2059.  
  2060. size (A);
  2061. size ([A]);
  2062. size (3.14);
  2063.  
  2064. size (Z);
  2065. size ([Z]);
  2066. size (3.14j);
  2067.  
  2068. printf ("\tsize");
  2069.  
  2070. sizeof (A);
  2071. sizeof ([A]);
  2072. sizeof (3.14);
  2073.  
  2074. sizeof (Z);
  2075. sizeof ([Z]);
  2076. sizeof (3.14j);
  2077.  
  2078. printf ("\tsizeof");
  2079. printf ("\n");
  2080.  
  2081. sort (A);
  2082. sort ([A]);
  2083. sort (3.14);
  2084.  
  2085. sort (Z);
  2086. sort ([Z]);
  2087. sort (3.14j);
  2088.  
  2089. printf ("\tsort");
  2090.  
  2091. sqrt (A);
  2092. sqrt ([A]);
  2093. sqrt (3.14);
  2094.  
  2095. sqrt (Z);
  2096. sqrt ([Z]);
  2097. sqrt (3.14j);
  2098.  
  2099. printf ("\tsqrt");
  2100.  
  2101. srand ();
  2102. srand ("clock");
  2103. srand (1);
  2104.  
  2105. printf ("\tsrand");
  2106.  
  2107. strsplt (S);
  2108. printf ("\tstrsplt");
  2109.  
  2110. strtod ("string");
  2111. strtod (S);
  2112. strtod ("1.2");
  2113. strtod (["1.2", "1e3"]);
  2114.  
  2115. printf ("\tstrtod");
  2116.  
  2117. sum (A);
  2118. sum ([A]);
  2119. sum (3.14);
  2120.  
  2121. sum (Z);
  2122. sum ([Z]);
  2123. sum (3.14j);
  2124.  
  2125. printf ("\tsum");
  2126.  
  2127. tan (A);
  2128. tan ([A]);
  2129. tan (3.14);
  2130.  
  2131. tan (Z);
  2132. tan ([Z]);
  2133. tan (3.14j);
  2134.  
  2135. printf ("\ttan");
  2136. printf ("\n");
  2137.  
  2138. tmpnam ();
  2139. printf ("\ttmpnam");
  2140.  
  2141. type (A);
  2142. type ([A]);
  2143. type (Z);
  2144. type ([Z]);
  2145. type (3.14);
  2146. type (S);
  2147. type ("string");
  2148.  
  2149. printf ("\ttype");
  2150.  
  2151. zeros (3,3);
  2152. zeros ([3], [3]);
  2153. zeros ([3,3]);
  2154.  
  2155. printf ("\tzeros");
  2156.  
  2157. printf ("\n");
  2158.  
  2159. srand(1);
  2160. rand("default");
  2161.  
  2162. //
  2163. //-------------------------- Test rand () -----------------------------
  2164. //
  2165. rand ("normal", 5, 1);
  2166. xrand = rand(4000,1);
  2167.  
  2168. mean = function(x)
  2169. {
  2170.   local(m);
  2171.  
  2172.   m = size (x)[1];
  2173.   if( m == 1 ) 
  2174.   { 
  2175.     m = size (x)[2];
  2176.   }
  2177.  
  2178.   return sum( x ) / m;
  2179. };
  2180.  
  2181. std = function(x)
  2182. {
  2183.   local(i, m, s);
  2184.  
  2185.   if(class(x) != "num") { error("std() requires NUMERICAL input"); }
  2186.  
  2187.   m = x.nr;
  2188.   if( m == 1 ) 
  2189.   { 
  2190.     return sqrt( sum( (x - mean(x)) .^ 2 ) / (x.nc - 1) );
  2191.   else
  2192.     for( i in 1:x.nc) {
  2193.       s[i] = sqrt( sum( (x[;i] - mean(x[;i])) .^ 2 ) / (x.nr - 1) );
  2194.     }
  2195.     return s;
  2196.   }
  2197. };
  2198.  
  2199. if (!(mean (xrand) > 4.9 && mean (xrand) < 5.1)) 
  2200.   { error ("error in random"); }
  2201. if (!(std (xrand) > 0.9 && std (xrand) < 1.1))
  2202.   { error ("error in random"); }
  2203. printf("\tpassed rand test...\n");
  2204.  
  2205. rand("default");
  2206.  
  2207. //
  2208. //-------------------------- Test norm () -----------------------------
  2209. //
  2210.  
  2211. tn = [1,2,3,4;2,1,2,3;3,2,1,2;4,3,2,1  ];
  2212. if (norm(tn,"m") != 4 ) { error ("incorrect norm computation"); }
  2213. if (norm(tn,"1") != 10) { error ("incorrect norm computation"); }
  2214. if (norm(tn,"i") != 10) { error ("incorrect norm computation"); }
  2215. printf("\tpassed norm test...\n");
  2216.  
  2217. //
  2218. //-------------------------- Test qr () ------------------------------
  2219. //
  2220.  
  2221. a = ohess(4);
  2222. qa = qr (a);
  2223. if (max (max (abs (qa.q*qa.r - a)))/(X*norm (a)*a.nr) > eps)
  2224.   { error ("possible qr() problems"); }
  2225.  
  2226. z = ohess (4) + ohess(4)*1i;
  2227. qz = qr (z);
  2228. if (max (max (abs (qz.q*qz.r -  z)))/(X*norm (z)*z.nr) > eps)
  2229.   { error ("possible qr() problems"); }
  2230.  
  2231. printf("\tpassed qr test...\n");
  2232.  
  2233. //
  2234. //-------------------------- Test schur () ----------------------------
  2235. //
  2236.  
  2237. a = randsvd (10, 10);
  2238. sa = schur (a);
  2239. if (max (max (abs (sa.z*sa.t*sa.z' - a)))/(X*norm (a)*a.nr) > eps)
  2240.   { error ("possible schur() problems"); }
  2241.  
  2242. z = rand (4,4) + rand(4,4)*1i;
  2243. sz = schur (z);
  2244. if (max (max (abs (sz.z*sz.t*sz.z' - z)))/(X*norm (z)*z.nr) > eps)
  2245.   { error ("possible schur() problems"); }
  2246.  
  2247. printf("\tpassed schur test...\n");
  2248.  
  2249. //
  2250. //-------------------------- Test schord () ----------------------------
  2251. //
  2252.  
  2253. s2a = schord (sa, 2, 4);
  2254. if (max (max (abs (s2a.z*s2a.t*s2a.z' - a)))/(X*norm (a)*a.nr) > eps)
  2255. { 
  2256.   error ("possible schord() problems"); 
  2257. }
  2258.  
  2259. s2z = schord (sz, 3, 1);
  2260. if (max (max (abs (s2z.z*s2z.t*s2z.z' - z)))/(X*norm (z)*z.nr) > eps)
  2261. { 
  2262.   error ("possible schord() problems"); 
  2263. }
  2264.  
  2265. printf("\tpassed schord test...\n");
  2266.  
  2267. //
  2268. //-------------------------- Test chol () -----------------------------
  2269. //
  2270.  
  2271. c = lehmer(10);
  2272. u = chol (c);
  2273. if (max (max (abs (u'*u - c)))/(X*norm (c)*c.nr) > eps)
  2274. { 
  2275.   error ("possible chol() problems"); 
  2276. }
  2277.  
  2278. cz = lehmer(10) + lehmer(10)*1j;
  2279. cz = symm(cz);
  2280. uz = chol (cz);
  2281. if (max (max (abs (uz'*uz - cz)))/(X*norm (cz)*cz.nr) > eps)
  2282. { 
  2283.   error ("possible chol() problems"); 
  2284. }
  2285.  
  2286. printf("\tpassed chol test...\n");
  2287.  
  2288. //
  2289. //--------------------------- Test inv () --------------------------------
  2290. //
  2291.  
  2292. a = randsvd(10,10);
  2293. b = ones(10,1);
  2294. x = inv(a) * b;
  2295. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2296. { 
  2297.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2298.   printf ("\tA*X - B:\n");
  2299.   abs (a*x - b)
  2300.   error ("possible inv() problems\n");
  2301. }
  2302.  
  2303. az = randsvd(10,10) + randsvd(10,10)*1j;
  2304. bz = rand(10,1) + rand(10,1)*1j;
  2305. xz = inv (az)*bz;
  2306. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2307. { 
  2308.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2309.   printf ("\tA*X - B:\n");
  2310.   abs (az*xz - bz)
  2311.   error ("possible inv() problems\n");
  2312. }
  2313.  
  2314. printf("\tpassed inv test...\n");
  2315.  
  2316. //
  2317. //-------------------------- Test solve () -----------------------------
  2318. //
  2319.  
  2320. //
  2321. // Real - General case
  2322. //
  2323.  
  2324. srand();
  2325. a = randsvd(10,10);
  2326. b = ones(10,1);
  2327. x = solve (a,b);
  2328. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2329. { 
  2330.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2331.   printf ("\tA*X - B:\n");
  2332.   abs (a*x - b)
  2333.   error ("possible solve() problems\n");
  2334. }
  2335.  
  2336. //
  2337. // Real - Symmetric case
  2338. //
  2339.  
  2340. s = symm (randsvd(10,10));
  2341. b = ones(10,1);
  2342. x = solve (s,b);
  2343. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2344. { 
  2345.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2346.   printf ("\tA*X - B:\n");
  2347.   abs (s*x - b)
  2348.   error ("possible solve() problems\n");
  2349. }
  2350.  
  2351. //
  2352. // Complex - General  case
  2353. //
  2354.  
  2355. az = randsvd(10,10) + randsvd(10,10)*1j;
  2356. bz = rand(10,1) + rand(10,1)*1j;
  2357. xz = solve (az,bz);
  2358. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2359. { 
  2360.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2361.   printf ("\tA*X - B:\n");
  2362.   abs (az*xz - bz)
  2363.   error ("possible solve() problems\n");
  2364. }
  2365.  
  2366. //
  2367. // Complex - Symmetric case
  2368. //
  2369.  
  2370. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2371. bz = rand(10,1) + rand(10,1)*1j;
  2372. xz = solve (sz,bz);
  2373. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2374. { 
  2375.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2376.   printf ("\tA*X - B:\n");
  2377.   abs (sz*xz - bz)
  2378.   error ("possible solve() problems\n");
  2379. }
  2380.  
  2381. printf("\tpassed solve test...\n");
  2382.  
  2383. //
  2384. //-------------------------- Test factor() / backsub() -----------------------------
  2385. //
  2386.  
  2387. //
  2388. // Real - General case
  2389. //
  2390.  
  2391. srand();
  2392. a = randsvd(10,10);
  2393. b = ones(10,1);
  2394. f = factor (a);
  2395. x = backsub (f,b);
  2396. if (max (max (abs(a*x - b)))/(X*norm (a)*a.nr) > eps)
  2397. { 
  2398.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2399.   printf ("\tA*X - B:\n");
  2400.   abs (a*x - b)
  2401.   error ("possible factor/backsub problems\n");
  2402. }
  2403.  
  2404. //
  2405. // Real - Symmetric case
  2406. //
  2407.  
  2408. s = symm (randsvd(10,10));
  2409. f = factor (s);
  2410. x = backsub (f,b);
  2411. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2412. { 
  2413.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2414.   printf ("\tA*X - B:\n");
  2415.   abs (s*x - b)
  2416.   error ("possible factor/backsub problems\n");
  2417. }
  2418.  
  2419. s = symm (randsvd(10,10));
  2420. f = factor (s, "s");
  2421. x = backsub (f,b);
  2422. if (max (max (abs(s*x - b)))/(X*norm (s)*s.nr) > eps)
  2423. { 
  2424.   printf ("\tThe condition # of s: %d\n", 1/rcond (s));
  2425.   printf ("\tA*X - B:\n");
  2426.   abs (s*x - b)
  2427.   error ("possible factor/backsub problems\n");
  2428. }
  2429.  
  2430. //
  2431. // Complex - General  case
  2432. //
  2433.  
  2434. az = randsvd(10,10) + randsvd(10,10)*1j;
  2435. bz = rand(10,1) + rand(10,1)*1j;
  2436. fz = factor (az);
  2437. xz = backsub (fz,bz);
  2438. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2439. { 
  2440.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2441.   printf ("\tA*X - B:\n");
  2442.   abs (az*xz - bz)
  2443.   error ("possible factor/backsub problems\n");
  2444. }
  2445.  
  2446. az = randsvd(10,10) + randsvd(10,10)*1j;
  2447. bz = rand(10,1) + rand(10,1)*1j;
  2448. fz = factor (az, "g");
  2449. xz = backsub (fz,bz);
  2450. if (max (max (abs (az*xz - bz)))/(X*norm (az)*az.nr) > eps)
  2451. { 
  2452.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2453.   printf ("\tA*X - B:\n");
  2454.   abs (az*xz - bz)
  2455.   error ("possible factor/backsub problems\n");
  2456. }
  2457.  
  2458. //
  2459. // Complex - Symmetric case
  2460. //
  2461.  
  2462. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2463. fz = factor(sz);
  2464. xz = backsub (fz,bz);
  2465. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2466. { 
  2467.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2468.   printf ("\tA*X - B:\n");
  2469.   abs (sz*xz - bz)
  2470.   error ("possible factor/backsub problems\n");
  2471. }
  2472.  
  2473. sz = symm (randsvd(10,10) + randsvd(10,10)*1j);
  2474. fz = factor(sz, "s");
  2475. xz = backsub (fz,bz);
  2476. if (max (max (abs (sz*xz - bz)))/(X*norm (sz)*sz.nr) > eps)
  2477. { 
  2478.   printf ("\tThe condition # of sz: %d\n", 1/rcond (sz));
  2479.   printf ("\tA*X - B:\n");
  2480.   abs (sz*xz - bz)
  2481.   error ("possible factor/backsub problems\n");
  2482. }
  2483.  
  2484. printf("\tpassed factor/backsub test...\n");
  2485.  
  2486. //
  2487. //-------------------------- Test svd ()   -----------------------------
  2488. //
  2489.  
  2490. a = randsvd(10,10);
  2491. s = svd (a);
  2492. if (max (max (abs (s.u*diag(s.sigma)*s.vt - a)))/(X*norm(a)*a.nr) > eps)
  2493. {
  2494.   error ("possible svd() problems"); 
  2495. }
  2496.  
  2497. z = randsvd(10,10) + rand(10,10)*1j;
  2498. sz = svd (z);
  2499. if (max (max (abs (sz.u*diag(sz.sigma)*sz.vt - z)))/(X*norm(z)*z.nr) > eps)
  2500. { 
  2501.   error ("possible svd() problems"); 
  2502. }
  2503.  
  2504. printf("\tpassed svd test...\n");
  2505.  
  2506. //
  2507. //------------------------------ Test lu() ---------------------------------
  2508. //
  2509.  
  2510. static (swap);
  2511.  
  2512. lu = function ( A )
  2513. {
  2514.   local (i, l, u, pvt, x)
  2515.  
  2516.   if (A.nr != A.nc) { error ("lu() requires square A"); }
  2517.  
  2518.   x = factor (A, "g");  // Do the factorization
  2519.  
  2520.   //
  2521.   // Now create l, u, and pvt from lu and pvt.
  2522.   //
  2523.  
  2524.   l = tril (x.lu, -1) + eye (size (x.lu));
  2525.   u = triu (x.lu);
  2526.   pvt = eye (size (x.lu));
  2527.  
  2528.   //
  2529.   // Now re-arange the columns of pvt
  2530.   //
  2531.  
  2532.   for (i in 1:max (size (x.lu)))
  2533.   {
  2534.     pvt = pvt[ ; swap (1:pvt.nc, i, x.pvt[i]) ];
  2535.   }
  2536.   return << l = l; u = u; pvt = pvt >>;
  2537. };
  2538.  
  2539. //
  2540. //  In vector V, swap elements I, J
  2541. //
  2542.  
  2543. swap = function ( V, I, J )
  2544. {
  2545.   local (v, tmp);
  2546.   v = V;
  2547.   tmp = v[I];
  2548.   v[I] = v[J];
  2549.   v[J] = tmp;
  2550.   return v;
  2551. };
  2552.  
  2553. a = randsvd(10,10);
  2554. lua = lu (a);
  2555. if (max (max (abs(a - lua.pvt*lua.l*lua.u)))/(X*norm (a)*a.nr) > eps)
  2556. { 
  2557.   printf ("\tThe condition # of a: %d\n", 1/rcond (a));
  2558.   printf ("\tA - p*l*u:\n");
  2559.   abs (a - lua.pvt*lua.l*lua.u)
  2560.   error ("possible lu()/factor() problems\n");
  2561. }
  2562.  
  2563. //
  2564. // Real
  2565. az = randsvd(10,10) + randsvd(10,10)*1j;
  2566. luz = lu (az);
  2567. if (max (max (abs (az - luz.pvt*luz.l*luz.u)))/(X*norm (az)*az.nr) > eps)
  2568. { 
  2569.   printf ("\tThe condition # of z: %d\n", 1/rcond (az));
  2570.   printf ("\tA - p*l*u:\n");
  2571.   abs (az - luz.ovt*luz.l*luz.u)
  2572.   error ("possible lu()/factor()() problems\n");
  2573. }
  2574.  
  2575. printf("\tpassed lu/factor test...\n");
  2576.  
  2577.  
  2578. //
  2579. //-------------------------- Test hess ()  -----------------------------
  2580. //
  2581.  
  2582. a = randsvd(10,10);
  2583. h = hess (a);
  2584. if (max (max (abs (h.p*h.h*h.p' - a)))/(X*norm (a)*a.nr) > eps)
  2585. { 
  2586.   error ("possible hess() problems");
  2587. }
  2588.  
  2589. z = randsvd(10,10) + randsvd(10,10)*1j;
  2590. hz = hess (z);
  2591. if (max (max (abs (hz.p*hz.h*hz.p' - z)))/(X*norm (z)*z.nr) > eps)
  2592. { 
  2593.   error ("possible hess() problems"); 
  2594. }
  2595.  
  2596. printf("\tpassed hess test...\n");
  2597.  
  2598. //
  2599. //-------------------------- Test lyap () ------------------------------
  2600. //
  2601.  
  2602. lyap = function ( A, B, C )
  2603. {
  2604.   local (X, sa, sb, tc)
  2605.  
  2606.   if (!exist (B)) 
  2607.   { 
  2608.     B = A'; // Solve the special form: A*X + X*A' = -C
  2609.   }
  2610.  
  2611.   if ((A.nr != A.nc) || (B.nr != B.nc) || (C.nr != A.nr) || (C.nc != B.nr)) {
  2612.     error ("Dimensions do not agree.");
  2613.   }
  2614.  
  2615.   //
  2616.   // Schur decomposition on A and B
  2617.   //
  2618.  
  2619.   sa = schur (A);
  2620.   sb = schur (B);
  2621.  
  2622.   //
  2623.   // transform C
  2624.   //
  2625.  
  2626.   tc = sa.z' * C * sb.z;
  2627.  
  2628.   X = sylv (sa.t, sb.t, tc);
  2629.  
  2630.   //
  2631.   // Undo the transformation
  2632.   //
  2633.  
  2634.   X = sa.z * X * sb.z';
  2635.  
  2636.   return X;
  2637. };
  2638.  
  2639. srand();
  2640. a = randsvd (10,10);
  2641. b = rand (10,10);
  2642. c = rand (10,10);
  2643.  
  2644. x = lyap (a, b, c);
  2645. if (max (max (abs (a*x + x*b + c)))/(X*norm(c)*norm(a)*norm(b)) > eps)
  2646. { 
  2647.   error ("possible problems with lyap() or sylv()"); 
  2648. }
  2649.  
  2650. printf("\tpassed lyap test...\n");
  2651.  
  2652. //
  2653. //-------------------------- Test eig () ------------------------------
  2654. //
  2655.  
  2656. trace = function(m) 
  2657. {
  2658.   local(i, tr);
  2659.  
  2660.   if(m.class != "num") { 
  2661.     error("must provide NUMERICAL input to trace()");
  2662.   }
  2663.  
  2664.   tr = 0;
  2665.   for(i in 1:min( [m.nr, m.nc] )) {
  2666.     tr = tr + m[i;i];
  2667.   }
  2668.  
  2669.   return tr;
  2670. };
  2671.  
  2672. eye = function( m , n ) 
  2673. {
  2674.   local(i, N, new);
  2675.  
  2676.   if (!exist (n))
  2677.   {
  2678.     if(m.n != 2) { error("only 2-el MATRIX allowed as eye() arg"); }
  2679.     new = zeros (m[1], m[2]);
  2680.     N = min ([m[1], m[2]]);
  2681.   else
  2682.     if (class (m) == "string" || class (n) == "string") {
  2683.       error ("eye(), string arguments not allowed");
  2684.     }
  2685.     if (max (size (m)) == 1 && max (size (n)) == 1)
  2686.     {
  2687.       new = zeros (m[1], n[1]);
  2688.       N = min ([m[1], n[1]]);
  2689.     else
  2690.       error ("matrix arguments to eye() must be 1x1");
  2691.     }
  2692.   }
  2693.   for(i in 1:N)
  2694.   {
  2695.     new[i;i] = 1.0;
  2696.   }
  2697.   return new;
  2698. };
  2699.  
  2700. //
  2701. // Standard eigenvalue problem
  2702. //
  2703.  
  2704. a = randsvd(10,10);
  2705. ta = trace (a);
  2706. sa = symm (a);
  2707. tsa = trace (sa);
  2708.  
  2709. z = randsvd(10,10) + randsvd(10,10)*1i;
  2710. tz = trace (z);
  2711. sz = symm (z);
  2712. tsz = trace (sz);
  2713.  
  2714. tol = 1.e-6;
  2715.  
  2716. if (!(ta < sum(eig(a).val) + tol && ta > sum(eig(a).val) - tol))
  2717. {
  2718.   error ("error in eig");
  2719. }
  2720. if (!(tsa < sum(eig(sa).val) + tol && tsa > sum(eig(sa).val) - tol))
  2721. {
  2722.   error ("error in eig");
  2723. }
  2724. if (abs(tz)+tol < abs(sum(eig(z).val)) && abs(tz)+tol > abs(sum(eig(z).val)))
  2725. {
  2726.   error ("error in eig");
  2727. }
  2728. if (abs(tsz)+tol < abs(sum(eig(sz).val)) && abs(tsz) > abs(sum(eig(sz).val)))
  2729. {
  2730.   error ("error in eig");
  2731. }
  2732.  
  2733. //
  2734. // Generalized eigenvalue problem
  2735. //
  2736.  
  2737. b = randsvd(10,10);
  2738. sb = symm (b) + eye(size(b))*3;
  2739. tb = trace (b);
  2740. tsb = trace (sb);
  2741.  
  2742. zb = randsvd(10,10) + randsvd(10,10)*1i;
  2743. szb = symm (zb) + eye(size(zb))*3;
  2744. tzb = trace (zb);
  2745. tszb = trace (szb);
  2746.  
  2747. eig(a,b);   // not sure of a good way to check these yet
  2748. eigs(sa,sb);
  2749. eigs(sa,sb);
  2750.  
  2751. eig(z, zb);
  2752. eigs(sz, szb);
  2753. eigs(sz, szb);
  2754.  
  2755. printf("\tpassed eig test...\n");
  2756.  
  2757. //
  2758. //-------------------------- Test fft () -----------------------------
  2759. //
  2760.  
  2761. if (100 != fft(ones(100,1))[1]) { error ("error in fft()"); }
  2762. printf("\tpassed fft test...\n");
  2763.  
  2764. //
  2765. //------------------------- Fibonacci Test -------------------------
  2766. //
  2767. //  Calculate Fibonacci numbers
  2768. //
  2769.  
  2770. i=1; 
  2771. while ( i < 2 ) { 
  2772.   i=i+1;
  2773.   a=0; b=1;
  2774.   while ( b < 10000 ) {
  2775.     c = b;
  2776.     b = a+b;
  2777.     a = c;
  2778.   }
  2779. }
  2780. if ( b != 10946 ) {
  2781.   error("failed fibonacci test");
  2782. else
  2783.   printf("\tpassed fibonacci test...\n");
  2784. }
  2785.  
  2786. //
  2787. //------------------------- Factorial Test -------------------------
  2788. //
  2789.  
  2790. fac = function(a) 
  2791. {
  2792.   if(a <= 1) 
  2793.   {
  2794.     return 1
  2795.   else
  2796.     return a*$self(a-1)
  2797.   }
  2798. };
  2799.  
  2800. if(fac(10) != 3628800)
  2801. { 
  2802.   error(); else printf("\tpassed factorial test...\n"); 
  2803. }
  2804.  
  2805. //
  2806. //--------------------------- ACK Test ----------------------------
  2807. //
  2808.  
  2809. ack = function(a, b) 
  2810. {
  2811.   if(a == 0) { return b + 1; }
  2812.   if(b == 0) { return $self(a-1, 1); }
  2813.   return $self(a-1, $self(a, b-1));
  2814. };
  2815.  
  2816. if(ack(2,2) != 7) 
  2817. {
  2818.   error("error in ack() test");
  2819.   else
  2820.   printf("\tpassed ACK test...\n");
  2821. }
  2822.  
  2823. //
  2824. //------------------------- Prime Test -----------------------------
  2825. //
  2826. // An example that finds all primes less than limit
  2827. //
  2828.  
  2829. primes = function (limit) 
  2830. {
  2831.   local(prime, cnt, i, j, k);
  2832.   
  2833.   i = 1; cnt = 0;
  2834.   for(k in 2:limit) 
  2835.   {
  2836.     j = 2;
  2837.     while(mod(k,j) != 0) 
  2838.     {
  2839.       j++;
  2840.     }
  2841.     if(j == k)             // Found prime
  2842.     {
  2843.       cnt++;
  2844.       prime[i;1] = k;
  2845.       i++;
  2846.     }
  2847.   }
  2848.   return prime;
  2849. };
  2850.  
  2851. if(max(size(primes(100))) != 25) 
  2852. {  
  2853.   error("error in prime test");
  2854.   else
  2855.   printf("\tpassed prime test...\n");
  2856. }
  2857.  
  2858. //
  2859. //--------------------------- Fib Min Test -----------------------------
  2860. //
  2861. //  fibmin() will minimize an arbitrary function 
  2862. //  in 1D using Fibonacci search
  2863.  
  2864. f065 = function(x)
  2865. {
  2866.   return (x - 0.65) * (x - 0.65);
  2867. };
  2868.  
  2869. fib = function(x)
  2870. {
  2871.   local (n, a, b);
  2872.   
  2873.   a = 1;
  2874.   b = 1;
  2875.   if (x >= 2)
  2876.   {
  2877.     n = x - 1;
  2878.     for (n in n:1:-1)
  2879.     {
  2880.       c = b;
  2881.       b = a + b;
  2882.       a = c;
  2883.       n = n - 1;
  2884.     }
  2885.   }
  2886.   return b;
  2887. };
  2888.  
  2889. //  Minimize a 1D function using Fibonacci search
  2890. //  f = function to minimize
  2891. //  xlo = lower bound
  2892. //  xhi = upper bound
  2893. //  n = number of iterations (the bigger the more accurate)
  2894.  
  2895. fibmin = function(f, xlo, xhi, n) 
  2896. {
  2897.   local(a, b, x, y, ex, ey, k, lo, hi);
  2898.   
  2899.   lo = xlo;
  2900.   hi = xhi;
  2901.   k = n;
  2902.   for (k in k:2:-1)
  2903.   {
  2904.     a = fib(k - 2) / fib(k);
  2905.     b = fib(k - 1) / fib(k);
  2906.     x = lo + (hi - lo) * a;
  2907.     y = lo + (hi - lo) * b;
  2908.     ex = f(x);
  2909.     ey = f(y);
  2910.     if (ex >= ey)
  2911.     {
  2912.       lo = x;
  2913.       else
  2914.       hi = y;
  2915.     }
  2916.     //  printf("%d: (%g %g) %g %g %g %g\n",  k, a, b, lo, hi, ex, ey);
  2917.   }
  2918.   return (lo + hi) / 2;
  2919. };
  2920.  
  2921. //
  2922. // Simple example using above function to mimize f065. Answer is 0.65
  2923. //
  2924.  
  2925. x = fibmin(f065, 0, 1, 30); // printf("f(%g)=%g\n", x, f065(x));
  2926. if (abs(x - 0.65) > 1e-6)
  2927. {
  2928.   printf("x = %f\n", x);
  2929.   error("failed fibmin test");
  2930. }
  2931.  
  2932. printf("\tpassed fibmin test...\n");
  2933.  
  2934. //
  2935. //--------------------- Nasty Function Test ------------------------
  2936. //
  2937.  
  2938. printf("\tStarting Nasty Function Test...");
  2939. printf("\tthis will take awhile\n");
  2940. check = function( a, b, c, d, e, f, g, h ) {
  2941.   if ( a+b+c+d == e+f+g+h && ...
  2942.        a^2+b^2+c^2+d^2 == e^2+f^2+g^2+h^2 && ...
  2943.        a^3+b^3+c^3+d^3 == e^3+f^3+g^3+h^3 ) {
  2944.     return 1;
  2945.   else
  2946.     return 0;
  2947.   }
  2948. };
  2949.  
  2950. for(a in 8:10) {
  2951.   for(b in 7:(a-1)) {
  2952.     for(c in 6:(b-1)) {
  2953.       for(d in 5:(c-1)) {
  2954.         for(e in 4:(d-1)) {
  2955.           for(f in 3:(e-1)) {
  2956.             for(g in 2:(f-1)) {
  2957.               for(h in 1:(g-1)) {                  
  2958.               if(check( a, b, c, d,  e, f, g, h ) || ...
  2959.                      check( a, e, c, d,  b, f, g, h ) || ...
  2960.                      check( a, f, c, d,  e, b, g, h ) || ...
  2961.                      check( a, g, c, d,  e, f, b, h ) || ...
  2962.                      check( a, h, c, d,  e, f, g, b ) || ...
  2963.                      check( a, b, e, d,  c, f, g, h ) || ...
  2964.                      check( a, b, f, d,  e, c, g, h ) || ...
  2965.                      check( a, b, g, d,  e, f, c, h ) || ...
  2966.                      check( a, b, h, d,  e, f, g, c ) || ...
  2967.                      check( a, b, c, e,  d, f, g, h ) || ...
  2968.                      check( a, b, c, f,  e, d, g, h ) || ...
  2969.                      check( a, b, c, g,  e, f, d, h ) || ...
  2970.                      check( a, b, c, h,  e, f, g, d ) || ...
  2971.                      check( a, e, f, d,  b, c, g, h ) || ...
  2972.                      check( a, e, g, d,  b, f, c, h ) || ...
  2973.                      check( a, e, h, d,  b, f, g, c ) || ...
  2974.                      check( a, f, g, d,  e, b, c, h ) || ...
  2975.                      check( a, f, h, d,  e, b, g, c ) || ...
  2976.                      check( a, g, h, d,  e, f, b, c ) || ...
  2977.                      check( a, b, e, f,  c, d, g, h ) || ...
  2978.                      check( a, b, e, g,  c, f, d, h ) || ...
  2979.                      check( a, b, e, h,  c, f, g, d ) || ...
  2980.                      check( a, b, f, g,  e, c, d, h ) || ...
  2981.                      check( a, b, f, h,  e, c, g, d ) || ...
  2982.                      check( a, b, g, h,  e, f, c, d ) || ...
  2983.                      check( a, e, f, g,  e, f, g, h ) || ...
  2984.                      check( a, e, f, h,  e, f, g, h ) || ...
  2985.                      check( a, e, g, h,  e, f, g, h ) || ...
  2986.                      check( a, f, g, h,  e, f, g, h ) ) { cnt++; }
  2987.               }
  2988.             }
  2989.           }
  2990.         }
  2991.       }
  2992.     }
  2993.   }
  2994. }
  2995.  
  2996. if(1) {  // figure out the value of cnt, and check!
  2997.   printf("\tpassed nasty function test...\n");
  2998. else
  2999.   error();
  3000. }
  3001.  
  3002. //
  3003. //------------------ Test More Advanced Functions --------------------
  3004. //
  3005.  
  3006. printf( "\tStarting the lqr/ode test..." );
  3007. printf( "\tthis will take awhile\n" );
  3008.  
  3009. lqr = function( a, b, q, r ) 
  3010. {
  3011.   local( k, s,...
  3012.          m, n, mb, nb, mq, nq,...
  3013.          e, v, d );
  3014.  
  3015.   m = size(a)[1]; n = size(a)[2];
  3016.   mb = size(b)[1]; nb = size(b)[2];
  3017.   mq = size(q)[1]; nq = size(q)[2];
  3018.   
  3019.   if ( m != mq || n != nq ) 
  3020.   {
  3021.     fprintf( "stderr", "A and Q must be the same size.\n" );
  3022.     quit
  3023.   }
  3024.     
  3025.   mr = size(r)[1]; nr = size(r)[2];
  3026.   if ( mr != nr || nb != mr ) 
  3027.   {
  3028.     fprintf( "stderr", "B and R must be consistent.\n" );
  3029.     quit
  3030.   }
  3031.     
  3032.   nn = zeros( m, nb );
  3033.     
  3034.   // Start eigenvector decomposition by finding eigenvectors of Hamiltonian:
  3035.     
  3036.   e = eig( [ a, solve(r',b')'*b'; q, -a' ] );
  3037.   v = e.vec; d = e.val;
  3038.     
  3039.   index = sort( real( d ) ).ind;
  3040.   d = real( d[ index ] );
  3041.   
  3042.   if ( !( d[n] < 0 && d[n+1] > 0 ) ) 
  3043.   {
  3044.     fprintf( "stderr", "Can't order eigenvalues.\n" );
  3045.     quit
  3046.   }
  3047.   
  3048.   chi = v[ 1:n; index[1:n] ];
  3049.   lambda = v[ (n+1):(2*n); index[1:n] ];
  3050.   s = -real(solve(chi',lambda')');
  3051.   k = solve( r, nn'+b'*s );
  3052.   
  3053.   return << k=k; s=s >>;
  3054.   
  3055. };
  3056.  
  3057. // Now run a little test problem.
  3058.  
  3059. k = 1; m = 1; c = .1;
  3060. a = [0     ,1    ,0    , 0;
  3061.     -k/m, -c/m,  k/m,  c/m;
  3062.      0,     0,    0,    1;
  3063.      k/m,  c/m, -k/m, -c/m ];
  3064. b = [ 0; 1/m; 0; 0 ];
  3065. qxx = diag( [0, 0, 100, 0] );
  3066. ruu = [1];
  3067. K = lqr( a, b, qxx, ruu ).k;
  3068.  
  3069. dot = function( t, x ) 
  3070. {
  3071.   global (a, b, K)
  3072.   return (a-b*K)*x + b*K*([1,0,1,0]');
  3073. };
  3074.  
  3075. x = ode ( dot, 0, 15, [0,0,0,0], .02, 1e-5, 1e-5 );
  3076.  
  3077. m = maxi( x[;2] );
  3078.  
  3079. if ( (abs( x[m;2] - 1.195 ) > 0.001)  || ...
  3080.      any (abs( x[x.nr;2,4] - 1 ) > 0.001) ) 
  3081. {
  3082.   printf( "\tfailed***\n" );
  3083.   else
  3084.   printf( "\tpassed the lqr/ode test...\n" );
  3085. }
  3086.  
  3087. printf("Elapsed time = %i seconds\n", toc() );
  3088. "FINISHED TESTS"
  3089.